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ABSTRACT 

We report the statistical properties of stars, brown dwarfs and multiple systems obtained from 
the largest radiation hydrodynamical simulation of star cluster formation to date that resolves 
masses down to the opacity limit for fragmentation (a few Jupiter masses). The initial condi- 
tions are identical to those of previous barotropic calculations published by Bate, but this time 
the calculation is performed using a realistic equation of state and radiation hydrodynamics. 
The calculation uses sink particles to model 183 stars and brown dwarfs, including 28 binaries 
and 12 higher-order multiple systems, the properties of which are compared the results from 
observational surveys. 

We find that the radiation hydrodynamical/sink particle simulation reproduces many ob- 
served stellar properties very well. In particular, whereas using a barotropic equation of state 
produces more brown dwarfs than stars, the inclusion of radiative feedback results in a stellar 
mass function and a ratio of brown dwarfs to stars in good agreement with observations of 
Galactic star-forming regions. In addition, many of the other statistical properties of the stars 
and brown dwarfs are in reasonable agreement with observations, including multiplicity as a 
function of primary mass, the frequency of very-low-mass binaries, and general trends for the 
mass ratio and separation distributions of binaries. We also examine the velocity dispersion 
of the stars, the distributions of disc truncation radii due to dynamical interactions, and copla- 
narity of orbits and sink particle spins in multiple systems. Overall, the calculation produces a 
cluster of stars whose statistical properties are difficult to distinguish from observed systems, 
implying that gravity, hydrodynamics, and radiative feedback are the primary ingredients for 
determining the origin of the statistical properties of low-mass stars. 
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1 INTRODUCTION 

Understanding the origin of the statistical properties of stellar sys- 
tems is the fundamental goal of a complete theory of star formation. 
Much attention has been paid to the origin of the stellar initial mass 
function (IMF), and there are many models that have been proposed 
for its origin (see the review oP Boimell, Larson & Zinnecker|200"7| l. 
However, the statistical properties of stellar systems include much 
more than just the IMF. A non-exhaustive list also includes the star 
formation rate and efficiency, the structure and kinematics of stellar 
groups and clusters, the properties of multiple stellar systems, jets, 
and protoplanetary discs, and the rotation rates of stars. A com- 
plete model must be able to explain the origin of all the statistical 
properties of stellar systems, and how these depend on variations 
in environment and initial conditions. While simplified analytic or 
semi-analytic models may be useful for understanding the role that 
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different processes play in the origin of some stellar properties, nu- 
merical simulations are almost certainly necessary to help us un- 
derstand the full complexity of the star formation process. 

To investigate the origin of the statistical properties of stel- 
lar systems through numerical simulations of star formation is dif- 
ficult because it is necessary to use sufficient resolution to en- 
sure that the processes involved are accurately modelled while si- 
multaneously producing a large number of stars from which sta- 
tistical properties can be derived. One approach is to perform a 
large number of hydrodynamical calculations of star formation 



in small molecular cloud cores (e.g. Delgado-Donate, Clarke & 



Batel[2004( Delgado-Donate et al.|[2004 



Goodwin, Whitworth & 



Ward-Thompson 2004a be, 2006, Stamat ellos, Rubber & Whit- 
worth 2007 , Stamatellos & Whitworth 2009| . Each calculation may 
only produce a few stars, but from the large ensemble of simula- 
tions the statistical properties can be studied. However, such calcu- 
lations begin with an arbitrary population of dense cores for their 
initial conditions, which may or may not be a good representation 
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of real dense cores. They also neglect evolution of the cores due to 
external processes such as growth of the cores by accretion, interac- 
tions with the turbulent environment in which they are embedded, 
and interactions between cores and protostellar systems which may 
be particularly important in dense star-forming regions. 

An alternative is to perform single large-scale hydrodynamical 
calculations of molecular clouds that each produce large numbers 
of stars. In these calculations, dense cores may form and evolve 
self-consistently from hydrodynamical flows on larger scales, and 
interactions between dense cores and protostellar systems can oc- 
cur naturally. Such calculations can be divided into two types: those 
that resolve small {%, 5 AU) scales to capture the opacity limit for 
fragmentation l |Low & Lynden-Bell||1976{ |Rees|1976[ >, and those 
that do not. 

Hydrodynamical calculations that do not resolve small scales 
miss some of the star formation and most binaries and discs. They 
usually seek to investigate the origin of the IMF and/or other large- 
scale properties such as the star formation rate. Early simulations 
of thi s type (|Klessen, Burkert & B ate"1998' 'Klessen & Burkert| 
lOOOllBonnell et al.|2001||Klessen|2 001 , Klessen & Burk ertpOOlf 



BonneU & Bate|2002[[Bonnell, Bate & Vine 2003 1 used isothermal 



equations of state and produced large numbers of stars, but used 
sink particles ( [Bate, Bormell & Price|1995] l with radii of hundreds 
of AU. More recent calculations have the investigated effects of 
additional physical processes on the origin of the IMF such as ra- 
diative transfer fOffner et al.'2009', 'Urban, Martel & Evans"2010) 
|Krumhoiz", Klein & McKee 2011 ), but these calculations also em- 
ploy sink particles with accretion radii of greater than 100 AU. 

The first hydrodynamical calculation to resolve the opacity 
limit for fragmentation and begin to probe the statistical properties 
of stars and brown dwarfs used a barotropic equation of state and 
sink particles with accretion radii of 5 AU, thus resolving many 
discs and binary and multiple systems ( [Bate, BonneU & Bromm| 
|2002a|b| [2003) . This calculation was followed by other similar 
calculations that probed the dependence of stellar properties on 
the mean thermal Jeans mass in the molecular cloud, the thermal 
behaviour of the gas, and the initial turbulent motions ( [Bate &| 
|BonneU|2003] |Bate|2005l |2009b^ . For example, these calculations 
showed that, when using a barotropic equation of state, the char- 
acteristic stellar mass depends primarily on the initial mean ther- 
mal Jeans mass in the cloud and not, for example, the initial tur- 
bulent power spectrum. These calculations were followed by those 
of other groups that also modelled the formation of stellar groups 
while simultaneously resolving discs and binaries ( [Li et al.|2004| 
lOffner, Klein & McKee|2008t . Most recently, calculations that re- 
solve these small scales and produce stellar groups have also begun 
to include the effects of additional physical processes such as mag- 
netic fields jPrice & Batel|2008t , radiative transfer (Bate 2009c I, 
and both of these combined (IPrice & Bate 20091. Using radiative 



transfer is found to dramatically decrease the amount of fragmen- 
tation, increase the characteristic stellar mass, decrease the propor- 
tion of brown dwarfs ( |Bate|2009c||Offner et al.|2009^ , and weaken 
the dependence of the characteristic mass of the IMF on the ini- 
tial Jeans mass (Bate 2009c). The latter effect may help to explain 
why the IMF is not observed to be strongly dependent on initial 
conditions, at least in our Galaxy ( [Bastian, Covey & Meyer|2010[ l. 
Stronger magnetic fields are found to decrease the star formation 
rate ( [Price & Bate|200"8l |2009l. However, in each of these calcu- 
lations, only a few dozen stars and brown dwarfs were produced, 
making it difficult to compare statistical properties with observa- 
tions in any detail. 

Currently, the only published hydrodynamical calculations 



that resolve the opacity limit for fragmentation and produce large 
numbers of stars, brown dwarfs (> 100) are those of |Batelp009al l. 
Two calculations were performed of SOO-M© molecular clouds, 
one using sink particle accretion radii of 5 AU, and an identical cal- 
culation using accretion radii of 0.5 AU that was not followed as far. 
The calculations used a barotropic equation of state to model the 
opacity limit for fragmentation. The former calculation produced 
more than 1250 stars and brown dwarfs, including well over 100 
multiple systems, that for the first time allowed a detailed compar- 
ison of a wide range of stellar properties with observations. It was 
found that many of the stellar properties were in good agreement 
with observed properties. For example, multiplicity was found to be 
a strongly increasing function of primary mass, the median separa- 
tion of multiple systems was found to decrease with primary mass, 
the mass ratios of very-low-mass (VLM) binaries (primary masses 
< 0.1 M0) were found to favour near equal-masses, and the rela- 
tive orbital orientations of triple systems were found to be some- 
what aligned. The good agreement with the observed properties 
of multiple stellar systems implies that such properties may orig- 
inate primarily from dissipative A'^-body dynamics, and that other 
physical processes such as radiative transfer and magnetic fields 
may play less of a role. However, the calculations produced far too 
many brown dwarfs relative to stars compared with a typical Galac- 
tic IMF. 

In this paper, we repeat the 5OO-M0 calculations of |Bate| 
(!2009a), but we use a realistic equation of state and include ra- 
diative transfer as implemented in the smaller calculations of |Bate| 
( |2009c| l and [Price & Bate] ( [2009[ l. Our aim is to investigate the 
effect of the realistic equation of state and radiative feedback on 
the star formation process in more detail than was possible with 
the earlier smaller calculations. In particular, we wish to determine 
whether the inclusion of radiative feedback can produce a more re- 
alistic IMF than that obtained by [Bate[ ( |2009a[ l, but retain the good 
agreement that was found when comparing the statistical properties 
of multiple stellar systems with observations. 



2 COMPUTATIONAL METHOD 

The calculations presented here were performed using a three- 
dimensional smoothed particle hydrodynamics (SPH) code based 
on the original version of ,BenZ|(1990; Benz et al.^19 90), but sub- 
stantially modified as descri bed in[Bate et al. ( 1995^, [W hitehouse^ 
Bate & Monaghan] pOOS] !, [Whitehouse & Bate[ ( [2006| ), |Price & 
Bate ( 2007 1, and paralleUsed using both OpenMP and MPI. 

Gravitational forces between particles and a particle's near- 
est neighbours are calculated using a binary tree. The smoothing 
lengths of particles are variable in time and space, set iteratively 
such that the smoothing length of each particle h = 1.2{m/ pY^^ 
where m and p are the SPH particle's mass and density, respectively 
(see [Price & Monaghan|2007| for further details). The SPH equa- 
tions are integrated using a second-order Runge-Kutta-Fehlberg in- 
tegrator with individual time steps for each particle ( [Bate et aT 
1995 I. To reduce numerical shear viscosity, we use the Morris & 



Monaghan] ( [ 1 997| > artificial viscosity with varying between 0.1 
and 1 while /3v = 2qv (see also [Price & Monaghan|2005) . 



2.1 Equation of state and radiative transfer 

The calculations presented in this paper were performed using radi- 
ation hydrodynamics with an ideal gas equation of state for the gas 
pressure p — pTgTZ/fi, where Tg is the gas temperature, p, is the 
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mean molecular weight of the gas, and TZ. is the gas constant. The 
thermal evolution takes into account the translational, rotational, 
and vibrational degrees of freedom of molecular hydrogen (assum- 
ing a 3: 1 mix of ortho- and para-hydrogen; see B oley et al. 2007j l. It 
also includes molecular hydrogen dissociation, and the ionisations 
of hydrogen and helium. The hydrogen and helium mass fractions 
are X = 0.70 and Y — 0.28, respectively. The contribution of 
metals to the equation of state is neglected. 

For this composition, the mean molecular weight of the gas 
is initially ^ = 2.38. The original calculation of |Bate|(2009a^ was 
performed using a barotropic equation of state which took the mean 
molecular weight of the gas to be /i — 2.46 and the initial temper- 
ature to be 10 K. To keep the same initial conditions (i.e. the same 
initial thermal energy of the gas), we set the initial temperature of 
the radiation hydrodynamical calculation to be 10.3 K. 

Two temperature (gas and radiation) radiative transfer in the 
flux-limited diffusion approximation is implemented as described 
by |Whitehouse et al. 1 2005 1 and Whi tehouse & Bate ( 2006 1, except 
that the standard explicit SPH contributions to the gas energy equa- 
tion due to the work and artificial viscosity are used when solving 
the (semi-)implicit energy equations to provide better energy con- 
servation. Energy is generated when work is done on the gas or ra- 
diation fields, radiation is transported via flux-limited diffusion and 
energy is transferred between the gas and radiation fields depend- 
ing on their relative temperatures, and the gas density and opacity. 
The gas and dust temperatures are assumed to be the same. Taking 
solar metallicity gas, the opacity is set to be the maximum of the in- 
terstellar grain opacity tables of |Pollack et al.|p985^ and, at higher 
temperatures when the dust has been destroyed, the gas opacity ta- 
bles of | Alexander |l |1975^ (the IVa King model) (see |Whitehouse &] 
|Bate|2006 l for further details). 

The cloud has a free boundary. To provide a boundary con- 
dition f or the radiative transfer we use the same method as IBatel 
1 2009c I. All particles with densities less than lO"'^^ g cm~'^ have 



their gas and radiation temperatures set to the initial values of 10.3 
K. This gas is two orders of magnitude less dense that the initial 
cloud (see Section [2!4| and, thus, these boundary particles surround 
the region of interest in which the star cluster forms. 

2.2 Sink particles 

Using the above realistic equation of state and radiation hydrody- 
namics means that as the gas collapses, each of the phases of pro- 
tostar formation are captured (Larson 1969). The initial collapse of 
a dense core proceeds almost isothermally, until the compressional 
rate heating of the gas exceeds the rate at which the gas can cool. 
At this point the collapse stalls, and a pressure supported fragment 
forms which [Larson |termed the 'first hydrostatic core'. The typical 
initial size and mass of this object is ~ 5 AU in radius and a few 
Jupiter masses (Mj). This first core accretes gas from the infalling 
envelope and if it is rotating rapidly it may undergo rotational dy- 
namical instabilities to form a disc ( |Bate|1998] l. Eventually, due to 
mass accretion (|Larson(1969f, dynamical instability ( |Bate|1998] ), 
and/or cooling (Tomida et al. 2010) the central temperature exceeds 
« 2000 K and molecular hydrogen begins to dissociate, absorbing 
thermal energy and resulting in a second collapse ( jLarson|[l969^ 
within the first core/disc. This collapse is halted when the dissoci- 
ation is complete and the 'second' or 'stellar core' forms ( |Larson| 
|1969[ l. This object initially has a radius of ~ 2 R© and a mass of 
~ 1 — 2 Mj. Subsequently it accretes to higher masses from the 
surrounding first core/disc and envelope. 

The calculations presented here include the physics necessary 



to follow each of these stages of protostar formation and evolu- 
tion. Indeed, the same code has been used to study the formation 
and evolution of first cores, pre-stellar discs, and stellar cores ( |Bate| 
|2010a[ |2011| ). However, as the collapse proceeds, the timesteps 
required to obey the Courant-Friedrichs-Levy (CFL) criterion be- 
come smaller and smaller Because we wish to evolve the large 
scales over timescales of ~ 10^ years, we cannot afford to follow 
the small scales (e.g. the stellar cores themselves). 

Instead, we follow the evolution of each protostar through the 
first core phase and into the second collapse (which beg ins at den- 
sities of 

3 



10 ^ g cm ''), but we insert a sink particle ([Bate et al. 



1995 I when the density exceeds 10 g cm . The timesteps re- 



quired to follow this evolution get very short, but the duration of 
the second collapse phase is quite brief and the use of individual 
particle timesteps mean that the calculation does not get slowed 
down for long. This density is just two orders of magnitude before 
the stellar core begins to form (density ~ 10"^ g cm^^), and a 
considerable improvement over previous similar barotropic calcu- 
lations. For example, [Bate et al.H2003^ and |Bate| ( [2009a^ inserted 
sink particles well before the onset of second collapse at densities 
of 10~" and 10"^° g cm , respectively. At these densities, sink 
particles might have been inserted before two fragments merged or 
before a fragment was disrupted. However, the time taken for a pro- 
tostar to evolve from 10"^ g cm~'^ to the formation of a stellar core 
is much less than a year (the free-fall time at this density is only one 
week!), so there is a no chance of protostellar fragments merging 
or begin disrupted between sink particle insertion and stellar core 
formation. 

In the main calculation discussed in this paper, a sink parti- 
cle is formed by replacing the SPH gas particles contained within 
'"acc = 0.5 AU of the densest gas particle in region undergoing sec- 
ond collapse by a point mass with the same mass and momentum. 
Any gas that later falls within this radius is accreted by the point 
mass if it is bound and its specific angular momentum is less than 
that required to form a circular orbit at radius race from the sink 
particle. Thus, gaseous discs around sink particles can only be re- 
solved if they have radii <; 1 AU. Sink particles interact with the 
gas only via gravity and accretion. There is no gravitational soften- 
ing between sink particles. The angular momentum accreted by a 
sink particle is recorded but plays no further role in the calculation. 

Since all sink particles are created within pressure-supported 
fragments, their initial masses are several Mj, as given by the opac- 
ity limit for fragmentation l |Low & Lynden-Bell|1976[|Rees|1976l ). 
Subsequently, they may accrete large amounts of material to be- 
come higher-mass brown dwarfs (^ 75 Mj) or stars (^ 75 Mj), but 
all the stars and brown dwarfs begin as these low-mass pressure- 
supported fragments. 

Sink particles are permitted to merge in the calculation if they 
passed within 0.01 AU of each other (i.e., ~ 2 R0). However, no 
mergers occurred during the calculation. 

2.3 Limitations of the sink particle approximation 

The benefits and potential problems associated with introducing 
sink particles into barotropic star formation calculations performed 
using SPH have been discussed by [Bate et aT] ( |1995^ , [Bate et aT] 
1I2OOT) and Bate'l i2009al l. Some of these problems are avoided in the 
calculation presented here. As mentioned above, it is no longer pos- 
sible that a fragment which has been replaced sink particle might 
have merged or been disrupted before stellar core formation if it 
had not been replaced by a sink particle. Another problem, found 
by [Batel ( |2009a) , was that the eccentricities of binary stellar systems 
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Figure 1. The dependence of the star formation rate on the sink particle accretion radius, race- For previously published barotropic and radiation hydrody- 
namical calculations of star formation in a turbulence 5O-M0 cloud we plot: the total stellar mass (i.e. the mass contained in sink particles) versus time (left 
panel), the number of stars and brown dwarfs (i.e. the number of sink particles) versus time (centre panel), and the number of stars and brown dwarfs versus 
the total stellar mass (right panel). The different line types are for: a barotropic calculation using race = 5 AU (dotted line; Bate et al. 2003); and radiation 
hydrodynamical calculations with race = 5 AU (short-dashed line; Bate 2009c), race = 0.5 AU (solid line; Bate 2009c), and race = 0.05 AU (long-dashed 
line; performed for this paper). Time is measured from the beginning of the calculations in terms of the free-fall time (1.90 X 10^ yr). It can be seen that 
the main change in the star formation is captured when changing from a barotropic equation of state to a radiation hydrodynamical calculation (both with 
race = 5 AU). Reducing the accretion radius used for the radiation hydrodynamical calculations has a smaller effect each time the accretion radius is reduced. 



with separations 1 — 20 AU were too high if sink particles with 
radii race = 5 AU were used due to the absence of small-scale 
dissipation, but that this was corrected if the radius was reduced 
to 0.5 AU. Therefore, for the calculation presented here we use 
Tacc ~ 0.5 AU. However, other problems with using sink particles 
remain. For example, using race ~ 0.5 AU, dissipative interactions 
between protostars on length-scales ^ 1 AU are still neglected. 
Similarly, discs smaller than ~ 1 AU in radius cannot be resolved. 

When using radiation hydrodynamics, there is a new problem 
— how to handle radiative feedback. If each stellar core itself was 



resolved (e.g Whitehouse & Bate 2006, Bate 2010a, 20111, the full 
radiative feedback from the protostars would be naturally modelled. 
However, introducing sink particles, means that the evolution inside 
the sink particle radius is neglected. In the simplest form where a 
sink particle consists only of a point mass and a hole, there is no 
radiative feedback from the stellar core and the inner part of an 
accretion disc on the rest of the calculation. This is the case for the 
calculations presented by |Bate||2009c[ l and also for the calculations 
presented in this paper. 

However, this means that the luminosity of the material inside 
the sink particle radius is omitted. There are three main sources of 
luminosity that may be inside the sink radius: the intrinsic luminos- 
ity of the protostar (or protostars), the luminosity generated by ac- 
cretion onto the protostar(s), and the luminosity of any other matter 
(e.g. an accretion disc). If the protostar accretes from a steady-state 
accretion disc, the accretion luminosity onto the protostar will dom- 
inate the luminosity of the disc itself. However, trying to determine 
the intrinsic luminosity of the protostar(s) and accretion luminosity 
is far from easy and we discuss the many problems below. 



_ GM,M, 
of a star of mass M* 



Hosokawa & Omukai 



(1) 

= 1 M (7) with a radius of i?, = 2 Rq (e.g. 
2009|) accreting at A/* = 1 x lO"*^ —-^ 



'i0yr 



is « 15 Lq whereas the luminosity of the stellar object itself is ~ 
I L© ( [Hosokawa & Omukai|2009) . As will be seen in Section [3.2.I| 
the typical accretion rates in the calculation presented here are an 
order of magnitude higher than this (~ 10""" Mq yr~^). There- 
fore, for the vast majority of the stars and brown dwarfs produced, 
the accretion luminosity dominates over stellar luminosity unless 
the accretion rate drops to very low levels. Very low accretion rates 
are usually obtained only after a star has had a dynamical inter- 
action and ejected from the dense star-forming cores, where upon 
its radiative feedback is no longer important for the subsequent star 
formation. Indeed, we define a star or brown dwarf whose accretion 
rate drops below 10~^ Mq yr~^ to have stopped accreting. 

For stars with masses ^ 3 Mq, whether accretion luminos- 
ity or ste llar luminosity dominates is more complex. Hosokawa &| 
Omukai |2009jl find that for accretion rates of 1 x lO"'' Mq yr~'^, 
accretion luminosity dominates for ^ 3 Mq, but with accretion 
rates of 1 x 10""^ Mq yr~^, accretion luminosity dominates for 
^ 9 Mq. There are only two stars with masses > 3 Mq pro- 
duced by the main calculation presented here: a 3.84 Mq star and 
a 3.38 Mq star. However, these two stars have average accretion 
rates > 5 x 10~® Mq yr~^, so the accretion luminosity is expected 
to dominate over stellar luminosity. 

We conclude that for low- and intermediate-mass star forma- 
tion like that considered in this paper, the intrinsic stellar luminosity 
can be confidently neglected. 



2.3.1 Intrinsic protostellar luminosity 

For a young low- or intermediate-mass (^ 3 Mq) protostar that is 
accreting at a significant rate (> 10^^ Mq yr^^), the accretion lu- 
minosity dominates and the intrinsic luminosity of the stellar object 
is negligible (e.g. |Offner et al.|2009[ [Hosokawa & Omukai|2009l l. 
For example, the accretion luminosity 



2.3.2 Accretion luminosity and sink particle accretion radii 

However, the accretion luminosity is a different issue. Because ac- 
cretion luminosity scales oc 1 / R (equation[TJ, excluding the radia- 
tion coming from inside a sink particle accretion radius means that 
the luminosity of an accreting stellar object is potentially underes- 
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timated by a factor of ~ R*/r^cc- Taking a protostellar radius of 
-R* « 3 — 5 Rq (typical for protos tars of mass 0.1 — 2 ac - 
creting at a rate of 10~^ yr~^; Hosokawa & Omukai|2009 1, 
this is a factor of ~ 200 — 300 when using an accretion radius of 
Tacc ~ 5 AU and a factor of ~ 20 — 30 when using race = 0.5 AU. 

Bate (2009c I performed radiation hydrodynamical calcula- 
tions of star formation that were similar to the calculation presented 
in this paper, but for smaller 5O-M0 clouds. He performed cal- 
culations using sink particle accretion radii of 5 AU and 0.5 AU 
(and, therefore, different fractions of the accretion luminosity) to 
determine the effect on the fragmentation. The initial conditions for 
these calculations were identical to the original barotropic calcula- 
tion of |Bate et al.| ( |2003} which used race = 5 AU. To investigate 
this issue further for this paper, we performed another calculation, 
identical to those of |Bate| ( [2009c| ), but using accretion radii of just 
0.05 AU (i.e. only « 2 — 3 times larger than a rapidly accreting low- 
mass star). In Fig.[T] we plot the total stellar mass (i.e. the mass in 
sink particles) and the number of sink particles as functions of time. 
We also plot the number of sink particles versus the total stellar 
mass. It is clear that the main effect of including radiative feedback 
is captured when going from a barotropic equation of state to radia- 
tion hydrodynamics and race ~ 5 AU which reduces the number of 
objects formed by more than a factor of 3. Decreasing the accretion 
radii from 5 AU to 0.5 AU to 0.05 AU has a progressively smaller 
and smaller effect. Although at t = 1.38 tg the calculation with 
T'acc = 0.05 AU had formed 11 objects while the calculation with 
face = 0.5 AU had formed 8 objects, the time evolution of the total 
stellar mass is almost identical for the two smallest accretion radii 
and until t = 1.35 tff (before the onset of the second burst of star 
formation; |Bate et al.|2003j l the number of stars only differs by one. 

Why should increasing the accretion luminosity by a factor of 
100 have such little effect compared to switching from a barotropic 
to radiation hydrodynamical calculation? There are three main rea- 
sons. First, the although the barotropic equation of state models the 
evolution of the maximum temperature/density during the collapse 
to form a protostar reasonably well, the temperature surrounding 
the object even before stellar core formation may be underestimated 
by up to an order of magnitude, particularly in a protostellar disc 
( [Whitehouse & Bate 2006 1. Thus, using radiation hydrodynamics 
rather than a barotropic equation of state dramatically decreases 
the propensity for massive discs to fragment even without stellar 
feedback. In |Bate et al.|p003| l, three-quarters of the brown dwarfs 
formed by disc fragmentation ( |Bate et al.|2002a t, so this has a par- 
ticularly important effect on the formation of low-mass objects. 
Second, in the envelope surrounding a protostar, the gas temper- 
ature depends on luminosity roughly as . Therefore, a large 
change in the luminosity does not actually translate into a large 
change in temperature. The Jeans mass scales with temperature as 
Tg'^^, so this translates into a change in the Jeans mass of ■ 
Therefore, changing the accretion radii from 5 AU, to 0.5 AU to 
0.05 AU, only decreases the Jeans mass near an existing protostar 
by factors of ~ 8, 3.4, and 1.4 compared to that expected from in- 
cluding the full accretion luminosity. Finally, protostars that form 
near each other usually form in a short period of time. Since when 
a new stellar core first forms most of the mass is still in the disc 
( |Bate|[l998l |2010a[ [Machida & Matsumoto|[20TT] >, the the accre- 
tion luminosity of the stellar core does not exceed the luminosity 
of the accreting first core/disc for some time. If the nearby frag- 
mentation occurs before the stellar luminosity becomes significant 
(e.g. |Bate|20l"T| >, the luminosity of the nearby protostar will be a 
good approximation despite the use of a sink particle. 



2.3.3 Problems with estimating sink particle luminosity 

Some studies have elected to try to include radiative feedback from 
inside the sink particle radius (e.g. |Offner et al.|2009[ [Urban et al.| 
[2010[[Krumholz et al.[201 1} . These investigations employ sink par- 
ticles with accretion radii in excess of 100 AU, so if nothing was 
done the radiative feedback from within this radius would be un- 
derestimated by a factor <; 10*. However, such attempts to include 
the missing radiative feedback rely on many approximations and 
assumptions. Problems include: deciding on the mass distribution 
within the accretion radius (for example, how much mass is in a 
disc versus a stellar object, or whether there is a single or multi- 
ple stellar system inside the accretion radius); deciding how much 
energy comes out in other forms of feedback (e.g. kinetic feed- 
back from jets, outflows, and winds); deciding how much energy 
is advected into the stellar object rather than radiated away; and 
deciding whether accretion is continuous or occurs in bursts. This 
list is not exhaustive, but it gives an indication of the difficulties 
associated with trying to accurately estimate the level of radiative 
feedback coming from inside the sink particle accretion radius. We 
now briefly discuss each of these problems and conclude that, gen- 
erally, the calculatio ns of[Offner et al.[p009) , p-ban et al.[ ( |20Tol l, 
and Knimh olz et al.[p01 1 1 likely overestimate the effects of radia- 
tive feedback. 

The simplest assumption is that all of the mass that enters a 
sink particle is immediately incorporated into a stellar object (e.g. 
Offner e^al. 2009; Urban et al. 2010; Krumholz et al. 201 1 1. How- 
ever, when using sink particles with sizes > 100 AU a consider- 
able fraction of the sink particle mass will still be in a protostellar 
disc. In fact, early in the star formation process, the vast major- 
ity of the mass will be in a disc rather than the stellar object. [Bate[ 
([1998)1 showed that the first cores formed from the collapse of ro- 
tating molecular cloud cores can in fact be pre-stellar discs (e.g. 
50 AU in radius) that may last thousands of years before a stel- 
lar core forms. Soon after stellar core formation, the disc mass can 
exceed the stellar mass by up to two orders of magnitude (Bate[ 
l998|[2010a[[Machida et al.|2010||Machida & Matsumoto 20lTy 



Bate|20lT) . Thus, allocating all of the mass of a sink particle to the 
stellar object can overestimate the level of radiative feedback soon 
after star formation by up to two orders of magnitude. 

Before stellar core formation, this overestimate can be even 
worse. During the thousands of years between first core/pre-stellar 
disc formation and the formation of a stellar core, the typical lumi- 
nosity of a first core/disc ranges from ~ 0.004 — 0.1 L0 depending 
on its rotation rate ( [Saigo & Tomisaka|2006[|20TT| . However, the 
mass within 100 AU can be substantial (e.g. up to 0.2 Mq; [Bate] 
[201 \) . Assuming that all of this mass is in a stellar core with a ra- 
dius « 3 R© and accreting at a rate of 1 x 10~^ solar masses per 
year (a typical accretion rate of a young object) gives an accretion 
luminosity of 20 L0 (i.e. 200 — 5000 times greater than the true 
value). This is particularly important because if fragmentation oc- 
curs near to an existing protostar it often occurs shortly after the 
first protostar formed. 

Later in the star formation process, it may perhaps be assumed 
that the mass of the disc is less than that of the stellar object (i.e. 
the stellar mass should be wrong by less than a factor of two; [Krat^ 
[ter & Matzner|2006[ [Kratter et al.|2008||20Tol l. However, even at 
this stage significant uncertainties remain. Another potentially big 
problem is that of multiple systems. For a low-mass sink particle, 
whether a sink particle contains a single object or a binary is not 
very important because, as mentioned above, the intrinsic stellar 
luminosity is usually negligible compared to the accretion luminos- 
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ity and, for a given accretion rate, all low-mass protostars are ex- 
pected to have similar radii. However, for more massive protostars 
the intrinsic stellar luminosity is expected to dominate over the ac- 
cretion luminosity. For zero-age-main-sequence stars with masses 
between 4 — 9 M0 (depending on accretion rate) and 20 M© the 
intrinsic stellar luminosity scales as oc . Furthermore, obser- 
vations show that massive stars preferentially have massive com- 
panions (e.g. jKobulnicky & Fryer|2007[ >. Therefore, if a high-mass 
sink particle actually contains an equal-mass binary rather than a 
single object, its luminosity will be overestimated by a factor of 
5.7 or, for an equal-mass triple system, a factor of 16. 
This problem does not occur in the simulations presented in this 
paper because they resolve the opacity limit for fragmentation and, 
thus, are expected to capture the formation of all stars and brown 
dwarfs. Furthermore, no massive protostars are formed. However, 
if sink particles with large accretion radii are used in calculations 
that do not resolve the opacity limit for fragmentation, then some 
massive sink particles may contain a multiple systems. 

Another problem with estimating the radiative feedback is de- 
termining the fraction of gravitational energy which is liberated as 
radiation rather than other forms of energy or retained within the 
accretion radius. In reality, some fraction of the energy will exit 
the sink particle accretion radius as kinetic energy of jets, outflows, 
and winds. [Offner & McKee|p011^ estimate that 1/4 of the gravi- 
tational energy may exit in kinetic form rather than as radiation for 
low-mass stars. For high-mass stars, the fraction is unknown. 

During the process of mass accretion onto a stellar object from 
a circumstellar disc, some of the energy will be advected into the 
star rather than all of it being radiated away. This efficiency fac- 
tor is not at all well understood. It is usually assumed that this 
fraction is very small (e.g. [Baraffe, Chabrier &~G allardo 20091, 
however if it is substantial this can lead to very different evolu- 
tionary paths of the stellar object (e.g . [Hartmann, Zhu & Calvet| 
[20TT] >. |Hosokawa, Offner & Krumholzl ( |201 1^ performed pre-mam- 
sequence stellar evolution calculations with either 'hot' or 'cold' 
accretion and found that the protostellar radius could be factors of 
2 — 4 larger with 'hot' accretion. This translates into an accretion 
luminosity that is 2 — 4 times lower than that predicted by 'cold' 
accretion models. 

The final problem we discuss here is one that was identi- 
fied more than two decades ago from observations of protostars. 
|Kenyon et al.|(T990| found that low-mass protostars are under lu- 
minous by two orders of magnitude when compared to simple es- 
timates of protostellar accretion. The survey of Eva ns et al.| j2009| 
recently confirmed this discrepancy, which is know as the "proto- 
stellar luminosity problem". Several solutions or partial solutions 
have been proposed to solve the problem (see [Offner & McKee] 
|201 1 ). As mentioned above, some of the gravitational potential en- 
ergy may be released in non-radiative forms (e.g. kinetic energy) 
and some may be advected into the stellar object leading to larger 
radii than usually assumed. Another possibility, first raised in the 
discovery paper ( [Kenyon et al.|I990| is that protostars do not ac- 
crete steadily, but that most of their mass is accreted in short bursts 
of accretion. During this time, the accretion luminosity would be 
very high, but most of the time they would be in a low luminosity 
state. Studies of such bursty accretion have become quite popu- 
lar recently (e.g. |Vorobyov & Basu' 2()05l|Zhu et al.|2009[|Baraffe| 
|et al., 2009 ). In terms of the effects of radiative feedback on star 
formation, if protostars spent the vast most of their time in a low- 
luminosity state, this would be similar to reducing the radiative 
feedback from the protostar to that emitted during the low luminos- 
ity state and ignoring the brief periods of high luminosity. Recently, 



this argument has been used by |Stamatellos, Whitworth & Hubber] 
([20TTJ to argue that those calculations that include continuous ra- 
diative feedback from sink particles may overestimate the effects 
of radiative feedback and, therefore, underestimate the amount of 
disc fragmentation and the formation of low-mass objects. 

2.3.4 Summary for sink particle luminosity 

In summary, there is no easy way to accurately include the radia- 
tive feedback from a sink particle on the rest of a radiation hydro- 
dynamical calculation. It is often assumed that including radiative 
feedback from inside the sink particle radius will be more accurate 
than excluding it ( [Offner et al.pOOgj [Urban et al.|2010[|Krumholz| 
[et al.|20fT| l. However, as the examples above show, this is far from 
certain. Indeed, as implemented in the literature, radiative feedback 
from sink particles almost certainly overestimates the effects of ra- 
diative feedback. Usually this overestimate is expected to be at the 
level of factors of ~ 2 — 4, but at the earliest stages of protostar 
formation the overestimate may be as much as 3 orders of magni- 
tude. This early radiative feedback would eventually be emitted by 
the source, but using a simple prescription it is emitted too early 
and may affect fragmentation at early times. Thus, all that can re- 
ally be said at the present time is that the actual affect of radiative 
feedback probably lies somewhere between the results obtained by 
excluding radiative feedback and the results that attempt to include 
radiative feedback. 

The choice made in this paper is to neglect radiative feedback 
coming from inside the sink particle radius, but to use as small 
an accretion radius as is computationally feasible. This model is 
elegant in that the results depend on a single parameter — the sink 
particle accretion radius. The effect of the missing radiation on the 
degree of fragmentation is tested by using smaller calculations with 
different accretion radii (Fig.[TJ. By decreasing the accretion radius 
by an order of magnitude (from 0.5 to 0.05 AU) we find that using 
sink particles with accretion radii of 0.5 AU may overestimate the 
degree of fragmentation by up to a factor of ~ 11/8 = 1.4. To put 
this in context, Knimholz et al. (201 1 1, who attempt to include sink 
particle luminosity, reduce their accretion radius by only a factor of 
two (from 400 AU to 200 AU) and produce 70% more stars. 

Naively, the calculations discussed in the remainder of this 
paper using sink particles with accretion radii of 0.5 AU underes- 
timate the accretion luminosity by factors of 20 — 30 (taking pro- 
tostellar accretion radii of 3 — 5 MQ; [Hosokawa & Omukai|2009| l. 
However, it has been pointed out (P. Andre, private communica- 
tion) that because protostars are observed to be under-luminous by 
1-2 orders of magnitude ( Kenyon et al.|1 990{ [Evans et al. 2009|, 
neglecting the radiative feedback from inside sink particles with ac- 
cretion radii of « 0.5 AU might be more accurate than including 
the 'missing' luminosity. If this is the case, then the level of radia- 
tive feedback in calculation discussed in the rest of this paper may 
be close to reality. 

2.4 Initial conditions 

The initial conditions for the calculation presented in this paper are 
identical to those of'Bate'f2009a'l. We followed the collapse of an 
initially uniform-density molecular cloud containing 500 Mq of 
molecular gas. The cloud's radius was 0.404 pc (83300 AU) giving 
an initial density of 1.2 x 1Q~^® g cm~^. At the initial temperature 
of 10.3 K, the mean thermal Jeans mass was 1 M© (i.e., the cloud 
contained 500 thermal Jeans masses). Although the cloud was uni- 
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Calculation Initial Gas Initial Accretion Gravity Time No. Stars No. Brown Mass of Stars & Mean Median 

Mass Radius Radii Softening Formed Dwarfs Formed Brown Dwarfs Mass Mass 

Mq pc AU AU tg Mq Mq Mq 



B2009a Main 


500 


0.404 


5 


4 


1.04 


^102 


^119 


32.6 


0.15 


0.06 












1.20 


> 235 


^ 355 


92.1 


0.16 


0.05 












1.50 


5S459 


^795 


191 


0.15 


0.06 


B2009a Rerun 


500 


0.404 


0.5 





1.04 


>94 


^164 


32.0 


0.12 


0.05 


Radiation Hydro 


500 


0.404 


0.5 





1.04 


> 50 


^ 19 


28.4 


0.41 


0.24 












1.20 


^ 147 


^ 36 


88.2 


0.48 


0.21 



Table 1. The parameters and overall statistical results for the two calculations of Bate (2009a) and the calculation presented here. The initial conditions for all 
calculations were identical. The Bate (2009a) calculations used a barotropic equation of state and the main calculation used sink particles with race = 5 AU 
and gravitational softening inside 4 AU, while the rerun calculation used race = 0.5 AU and no gravitational softening. The calculation presented here was 
identical to the rerun calculation, except that it used a realistic equation of state with radiation hydrodynamics. The calculations were run for different numbers 
of initial cloud free-fall times due to computational limitations. Brown dwarfs are defined as having final masses less than 0.075 Mq . The numbers of stars 
(brown dwarfs) aie lower (upper) limits because some brown dwarfs were still accreting when the calculations were stopped. Including radiative feedback 
decreases the number of objects formed at the same time by a factor of 3.2 — 3.7, and increases the mean and median masses of the objects by factors of 
3 and 4, respectively. The amount of gas converted into stars only decreases by 4 — 11% compared to the barotropic calculations at the same times. 



form in density, we imposed an initial supersonic 'turbulent' veloc- 
ity field in the same manner as |Ostriker, Stone & GamimeljlOOl) 
and [Bate et aL] (j2003j. We generated a divergence-free random 
Gaussian velocity field with a power spectrum P{k) oc k"*', where 
k is the wavenumber. In three dimensions, this results in a ve- 
locity dispersion that varies with distance. A, as u{\) oc A^^^ in 
agreement with the observed Larson scaling relations for molec- 
ular clouds (f Larson|1981| ). The velocity field was generated on a 
128^ uniform grid and the velocities of the particles were interpo- 
lated from the grid. The velocity field was normalised so that the 
kinetic energy of the turbulence was equal to the magnitude of the 
gravitational potential energy of the cloud. Thus, the initial root- 
mean-square (rms) Mach number of the turbulence was = 13.7. 
The initial free-fall time of the cloud was tg — 6.0 x 10^^ s or 
1.90 X 10^ years. 

Since the initial conditions for the calculation are identical to 
those of |Bate|(2009a[ l and including radiative transfer does not al- 
ter the temperature of the gas significantly until just before the 
first protostar forms, the early evolution of the cloud was not re- 
peated. Instead, the radiation hydrodynamical calculation was be- 
gun from a dump file from the original calculation just before the 
maximum density exceeded 10~^^ g cm~'^. This method was also 
used for the radiation hydrodynamical calculations presented by 
|Bate| ( (2009c^ , which were radiation hydrodynamical versions of the 
earlier barotropic calculations published by |Bate et al.| ( |2003] l and 
[Bate & BonneiilpOOS} . 



Supercomputer, an SGI Altix ICE 8200. Using 256 compute cores, 
it required 6.3 million CPU hours (i.e. 34 months of wall time). 



3 RESULTS 

The calculation presented here is a radiation hydrodynamical ver- 
sion of the barotropic calculations presented by Bate'(2009a[>. [Bate| 
presented the results from two calculations that differed only in 
the value of the sink particle accretion radius used and the gravita- 
tional softening between sink particles. The main calculation with 
face = 5 AU produced 1254 stars and brown dwarfs in 1.50tfi 
(285,350 years) and the rerun calculation used race = 0.5 AU 
and produced 258 objects in LOSStg (197,460 years). See Table 
[T]for a summary of the statistics, including the numbers of stars 
and brown dwarfs produced by the two calculations, the total mass 
that had been converted to stars and brown dwarfs, and the mean 
and median stellar masses. Bate ( 2009a j compared a large number 
of statistical properties of the stars and brown dwarfs formed in 
the calculations with observations, finding that many were in good 
agreement with observations (see the Introduction). In this paper, 
we use the same analysis methods as those used by |Bate| l |2009a[ l, 
and we compare and contrast the results both with the results from 
the earlier barotropic calculations and with the results of observa- 
tional surveys. 



3.1 The evolution of the star-forming cloud 



2.5 Resolution 

The local Jeans mass must be resolved throughout the calculation 
to model fragmentation corre ctly ([Ba te & Bu rkert|1997 Truelove 
let al. ^1997 Whitwo rth|1998l|Bosset al.|200dnH'ubber, Goodwin 



& Whitworth 2006). The original barotropic calculation of IBate 



l |2009a| used 3.5 x 10^ particles to model the 500-Mq cloud and 
resolve the Jeans mass down to its minimum value of ~ 0.0011 
(1.1 Mj) at the maximum density during the isothermal phase of the 
collapse, pcrit = 10~^^ g cm""^. Using radiation hydrodynamics 
results in temperatures at a given density no less than those given 
by the original barotropic equation of state (e.g. |Whitehouse & Bate| 
12006) and, thus, the Jeans mass is also resolved in the calculations 



presented here. 

The calculation was performed on the University of Exeter 



As mentioned in Section [24] the radiation hydrodynamical calcu- 
lation was not re-run from the initial conditions, but was started 
from the last dump file of the original barotropic calculation before 
the density exceeded 10~^^ g cm""*. Before this point the initial 
'turbulent' velocity field had generated density structure in the gas, 
some of which was collected into dense cores which had begun to 
collapse. Those readers interested in this early phase should refer 
to [Batel ( |2009a) for fi gures and further details. 

The radiation hydrodynamical calculation was started from 
t = 0.700 tff. In the barotropic calculation, the first sink parti- 
cle was inserted ai. t — 0.715 ts, some 2.9 x 10'^ years later. 
Using radiation hydrodynamics, the first sink particle is inserted at 
t = 0.727 tff . The slightly later time (2.2 x 10^ yrs) is primarily be- 
cause in the original calculation sink particles were inserted when 
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Figure 2. The global evolution of column density in the radiation hydrodynamical calculation. Shocks lead to the dissipation of the turbulent energy that 
initially supports the cloud, allowing parts of the cloud to collapse. Star formation begins at t = 0.727 tg in a collapsing dense core. By t = l.lOtff the 
cloud has produced six main sub-clusters, and by the end of the calculation these sub-clusters started to merge and dissolve. Each panel is 0.6 pc (123,500 AU) 
across. Time is given in units of the initial free-fall time, tff = 1.90 X 10^ yr. The panels show the logarithm of column density, A'^, through the cloud, with 
the scale covering— 1.4 < logA^ < 1.0 with Af measured in g cm~^. White dots represent the stars and brown dwarfs. 



Stellar, brown dwarf, and multiple star properties 9 








Figure 3. The global evolution of gas temperature in the radiation hydrodynamical calculation. Initially, the gas is almost isothermal on large-scales, but as 
groups of stars begin to form they locally heat the gas. Soon after t = 1.15%, the merger of two stellar groups at the centre of the cluster and increased 
accretion rates onto the stars heats the inner 0.2 pc of the cluster. Each panel is 0.6 pc (123,500 AU) across. Time is given in units of the initial free-fall time, 
tg = 1.90 X 10^ yr. The panels show the logarithm of mass weighted gas temperature, Tg, through the cloud, with the scale covering 9 — 50 K. White dots 
represent the stars and brown dwarfs. 



the density exceeded 10"^" g cm~^ (when the fragment was in the 
'first hydrostatic core' stage of evolution) whereas with the radi- 
ation hydrodynamics we do not insert sink particles until halfway 
through the second collapse phase at a density of 10~^ g cm^"' 
(see Section [2^ for further details). The first core phase typically 
lasts several thousand years, depending on the amount of rotation 
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In the panels of Figs. [2] and [5] we present snapshots of the 
global evolution of the calculation. Fig. |2] displays the column 
density (using a red-yellow- white colour scale), while Fig. [3] dis- 



plays the mass-weighted temperature in the cloud (using a blue-red- 
yellow-white colour scale). The column-density snapshots may be 
compared with the figures in |Bate| ( (2009a^ that show the barotropic 
evolution. Animations of both the main barotropic calculation and 
the radiation hydrodynamical calculation can be downloaded from 
|http://w WW. astro .ex . ac .uk/people/mbate/ ^ | 

|BateY2009c^ found that barotropic and radiation hydrodynam- 
ical calculations diverge quickly on small scales. In particular, it 
was found that many massive circumstellar discs that quickly frag- 
mented in barotropic calculations did not fragment, or took much 
longer to fragment, when using radiation hydrodynamics. The ac- 
cretion luminosity released as gas falls onto a disc and then spirals 
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Figure 4. The star formation rate obtained using a barotropic equation of state or radiation liydrodynamics for the SOO-Mq turbulent molecular cloud. We 
plot: the total stellar mass (i.e. the mass contained in sink particles) versus time (left panel), the number of stars and brown dwarfs (i.e. the number of sink 
particles) versus time (centre panel), and the number of stars and brown dwarfs versus the total stellar mass (right panel). The different line types are for: the 
main barotropic calculation using r'acc = 5 AU (dotted line; Bate 2009a); the rerun barotropic calculation using race = 0.5 AU (dashed line; Bate 2009a); 
and the radiation hydrodynamical calculation with r^cc = 0.5 AU (solid line). Time is measured from the beginning of the calculation in terms of the free-fall 
time of the initial cloud (bottom) or years (top). The rate at which mass is converted into stars is almost unaltered by the introduction of radiative feedback, 
but the number of stars and brown dwarfs is decreased by a factor of three. 



in towards the central protostar is often sufficient to heat a disc 
and prevent its fragmentation. An individual low-mass protostar 
can also produce substantial heating of the surrounding cloud out 
to thousands of AU (depending on the protostar's accretion rate) 
which can inhibit the fragmentation of nearby collapsing filaments. 
|Bate| ( [2009c| l found that the temperature field around rapidly accret- 
ing protostars could vary significantly on timescales of hundreds 
to thousands of years due to variations in the accretion rates of 
the protostars and their discs, particularly when protostars undergo 
dynamical interactions that perturb their discs. The same effects 
are found in the calculation presented here, though with an order 
of magnitude more objects it is impossible for us to compare and 
contrast the barotropic and radiation hydrodynamical evolutions of 
individual objects as was done by |BateH2009c| l. The temperature 
variability is best seen in an animation. 

The large-scale influence of the radiative feedback from the 
young protostars can be seen in Fig. [3] Each panel in this figure 
measure 0.6 pc (123,500 AU) across. As in the calculations o f|Bon-| 
|nell et al.||2003| > and |Bate|j2009a| l, the star formation in the cloud 
occurs in small groups, often formed within larger-scale filaments. 
Initially, each group contains only a few low-mass objects and the 
heating of the surrounding gas is limited to their immediate vicin- 
ity (a few thousand AU). However, as the stellar groups grow in 
number and some of the stars grow to greater masses, the heating 
can be seen to extend to larger and larger scales. At t ~ 1.15 tg, 
the merger of several stellar groups occurs near the centre of the 
cloud and the protostellar accretion rates also increase. This pro- 
duces a burst of radiation that heats the centre of the cloud out to 
distances of ~ 0.2 pc (~ 80, 000 AU). Several bursts between this 
time and the end of the calculation {t = 1.20 iff ) continuously heat 
the centre of the cloud. 

[Bate] ( |2009a[ > followed the main barotropic calculation to 
1.50 iff (285 350 yr) at which time 38% of the gas had been con- 
verted into 1254 stars and brown dwarfs. Unfortunately, due to the 
extra computational expense of resolving the gas near sink parti- 
cles to 0.5 AU and the implicit radiative transfer we are only able 
to follow the radiation hydrodynamical case to 1.20 tfi (228 280 
yr) which is 9.0 x lO** years after the star formation began. At this 



stage 88.2 M0 of gas (18%) has been converted into 183 stars and 
brown dwarfs. Table[T|gives the numbers of stars and brown dwarfs 
and their mean and median masses for the radiation hydrodynam- 
ical calculation, and for the barotropic calculations. The informa- 
tion is given at the end points of each of the calculations and, where 
possible, at the same times to allow direct comparison between the 
different calculations. 

As was found by |Bate| ( |2009c) for smaller 5O-M0 clouds (see 
also Fig. [TJ, the radiative feedback dramatically reduces the num- 
ber of objects formed. Comparing the main barotropic calculation 
and the radiation hydrodynamical calculation &Xt= 1.20 tg, the 
former had produced 590 objects while the latter has only produced 
183 (less than 1/3). However, this reduction in the number of ob- 
jects is not the same for all stellar masses. The main barotropic 
calculation produced 235 stars and 355 brown dwarfs in the same 
time that the radiation hydrodynamical calculation produced 147 
stars and 36 brown dwarfs. Thus, the number of stars has only been 
reduced to 63% of the barotropic value, but the number of brown 
dwarfs has been slashed to just 10%! This change in the stellar mass 
distribution is also reflected in the mean and median masses (Table 
^ with the mean mass increasing by a factor of three from 0. 16 M0 
to 0.48 M0 and the median mass increasing by a factor of 4.2 from 
0.05 M© to 0.21 M© (measured at 1.20 ts). The maximum stel- 
lar mass is 3.84 M0, whereas the main barotropic calculation had 
produced a star of mass 3.13 M© at the same time and went on to 
produce a star of 5.4 M© by the end of the calculation. We inves- 
tigate the change in the distribution of stellar masses further in the 
next section. Before that, in Fig. |4]we examine the star formation 
rate in terms of mass and number of stars and brown dwarfs. In the 
left panel, we plot the total stellar mass as a function of time for the 
barotropic calculations of |Bate| ( [2009a| and the radiation hydrody- 
namical calculation. It can be seen that in terms of stellar mass, 
there is a slow star formation rate of « 5 x 10~* Mq yr~^ from 
^ 0.8 — 1.0 tg followed by an increase to ^ 2 x 10"^ Mq yr~^ 
after ~ 1.0 iff . The star formation rate is quite constant after this 
transition and, despite the dramatic effect of the radiative feedback 
in heating the cloud (Fig. [3]l, there is no evidence of a decreasing 
rate near the end of the radiation hydrodynamical calculation. In 
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the main barotropic calculation, there is a hint of a decrease after 
1.40 iff. This is not surprising since after this point more than a 
third of the gas has been converted to stars and some of the remain- 
ing gas is drifting off to large distances. 

Comparing the barotropic and radiation hydrodynamical cal- 
culations, the star formation rate in terms of Mq yr^^ is almost 
entirely unaffected by the inclusion of radiative feedback. At the 
end of the radiation hydrodynamical calculation, 88.2 M© of gas 
has been converted to stars while, at the same time, 92.1 of 
gas had been converted to stars in the main barotropic calculation 
(a difference of only 4%). |Bate| P009c^ also found that radiative 
feedback had little effect on the rate at which gas was converted to 
stars — with one type of initial condition the rate decreased by 4%, 
while in the other it increased by 15%. Similarly, Krumholz et al.| 
(j20TTJ recently modelled star formation in a 1000 M© cloud and 
found very similar star formation rates in terms of Mq yr~^ with 
and without radiative feedback. 

We note that a general problem with hydrodynamical models 
of star formation in bound molecular clouds (whether they include 
radiative feedback or not) is that the star formation rate is much 
quicker than observed. The observed star formation efficiency per 
free-fall time is 3 — 6% ( [Evans et al.|2009| ), whereas the peak rate 
in the calculation presented here is 76% (i.e. an order of magnitude 
greater). The solution(s) to this problem may be that star forma- 
tion occurs in globally unbound molecular clouds ( Clark & Bonnell 
2004; Cl ark et ari [2005 ), or that magnetic support (Price & Bate 
2008, 2009) , kinetic feedback ( ,Matzner & McKee,2000„Krumholz, 
et al.| 2006[ |Nakamura & Li|[2007^ or a combination (e.g. |Wang 
et al.|2010 l reduce the star formation rate. Investigating these ef- 
fects is beyond the scope of this paper, but they certainly warrant 
future investigation. 

Rather than altering the rate at which gas is converted into 
stars, the effect of radiative feedback is to convert mass into fewer 
stars and brown dwarfs by inhibiting fragmentation of the gas. The 
reduction in the rate of production of new protostars is clear from 
the centre and right panels of Fig.|4] Throughout the evolution, the 
radiation hydrodynamical calculation consistently produces new 
objects at about 1/3 the rate of the barotropic calculations. How- 
ever, as with the rate at which mass is converted into stars, there 
is no evidence that the radiative heating of the central regions of 
the cloud is reducing the rate at which new stars are being formed. 
This is in contrast to the results obtained by |Krumholz et al.1 ( |201 
who found that the initial rate of protostar formation was similar 
with and without radiative feedback, but that as the calculation pro- 
gressed the rate of protostar formation dropped off much faster with 
radiative feedback than without. Part of this difference between the 
results here and those of |Krumholz et al.] is certainly due to the dif- 
ferent initial conditions. The initial conditions of IKrumholz et all 
are more than seven times denser (twice the mass within a cloud ra- 
dius of 0.26 pc rather than 0.4 pc) and they are centrally condensed 
rather than uniform in density. This produces qualitatively differ- 
ent protostar formation even without radiative feedback: Krumholz' 
|et al. find a monotonically decreasing protostar formation rate as a 
function of mass without feedback, whereas we obtain a constant 
rate (dotted line, right panel of Fig. |4]l. The centrally-condensed 
initial conditions of Kru mholz et al.l also favour the formation of 
massive stars ( jGirichidis et al.|2011[ l, while uniform initial condi- 
tions result in the formation of many sub-clusters ^Bonnell et al.l 
|2003[ l which only later merge into a single cluster (Bate 2009al. 
The higher density cloud of [Krumholz et al.] gives a star formation 
rate of 2.4 x 10"'^ yr~ , more than an order of magnitude 
greater than obtained here. The combination of a higher star forma- 
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Figure 5. The cumulative IMF produced at the end (t = 1.20 %) of the 
radiation hydrodynamical calculation (solid line), compared to the IMF ob- 
tained from the main barotropic calculation by Bate (2009a) at t = 1.20 tfi 
(long-dashed line) and t = 1.50 % (dotted line), and the rerun barotropic 
calculation at t = 1.04 iff (short-dashed line). We also plot using the 
dot-long-dashed line, the cumulative IMF from the parameterisation of 
the observed IMF by Chabrier (2005). The vertical dashed line marks the 
stellar/brown dwaii' boundary. The form of the barotropic IMF does not 
change substantially from 1.04 to 1.50 ff (Kolmogorov-Smimov tests show 
the three distributions to be indistinguishable), but introducing radiative 
feedback substantially increases the median stellar mass and changes the 
IMF from being dominated by brown dwarfs to being dominated by stars. 
A Kolmogorov-Smimov test gives less than 1 chance in 10^^ that the 
barotropic and radiation hydrodynamic IMFs at i = 1.20 iff are drawn 
from the same underlying distribution. On the other hand, a Kolmogorov- 
Smimov test gives a 50.5 percent probability that the radiation hydrody- 
namic IMF could have been drawn from Chabrier's fit to the observed IMF 
(i.e. the two mass functions are indistinguishable). 

tion rate and a bias towards massive star formation mean that when 
radiative feedback is included it has a much greater effect. Finally, 
as discussed in Section [23||Krumholz et alTj include radiative feed- 
back from sink particles which may be over-estimated in contrast 
to our feedback which may be under-estimated. 

If the radiation hydrodynamical calculation presented here 
was followed further then, as shown by |Bate| ( [2009a| l, the even- 
tual result would be that most of the small groups and sub-clusters 
would merge into a single cluster surrounded by a halo of ejected 
objects. Unfortunately, the calculation cannot be followed that 
far due to computational limitations. However, it is already clear 
by comparing the main barotropic calculation with the radiation 
hydrodynamical calculation at the same time, that the number 
of ejected objects is substantially lower with radiative feedback. 
This is because of the reduction in small-scale fragmentation. The 
ejected objects in the barotropic calculations primarily come from 
small dense multiple systems. With radiative heating, each multi- 
ple system or stellar group contains fewer objects so their are fewer 
dynamical interactions and ejections. 

3.2 The initial mass function 

In past barotropic calculations of the formation of stellar groups 
and clusters, the number of brown dwarfs produced has consistently 
been greater than the number of stars, in disagreement with obser- 
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Figure 6. The cumulative stellar mass distributions produced at various 
times during the radiation hydrodynamical calculation. The times are t = 
0.90 tff (dotted line), t = 1.00 tfl (short-dashed line), t = 1.10 % (long- 
dashed line), and t = 1.20 iff (solid line). The vertical dashed line marks 
the stellar/brown dwarf boundary. The form of the stellar mass distribution 
does not change significantly during the radiation hydrodynamical calcula- 
tion, though as more stars are formed the maximum stellar mass increases. 



vations of Galactic star-forming regions Pate et al.|2003l|Bate &| 
|BonnelI|2005| |Bate'2005! [2009a|b> . The radiation hydrodynamical 
calculations of Bate (2009cjl showed that radiative feedback pro- 
vides a potential solution to this over production of low-mass ob- 
jects. Although the clouds studied were an order of magnitude less 
massive than the calculation presented here, |Bate| |2009c| found 
that the number of objects formed when including radiative feed- 
back was reduced by a factor of ~ 3.8 compared to the barotropic 
calculation of [Bate et al.|P003] l and the mean and median masses 
of the objects increased by factors of ~ 4. The effects of radiative 
feedback found here are very similar: a reduction in the number of 
objects by a factor of 590/183 = 3.2 and increases of the mean 
and median masses by factors of ~ 3 and ~ 4, respectively (Table 
[T](. However, with an order of magnitude more objects we are able 
to examine the change in the mass function in more detail. 

In Fig. |5] we compare the cumulative IMF at the end of the 
radiation hydrodynamical calculation (solid line) with the IMFs of 
the main and rerun barotropic calculations ^Bate,2009a^ . With ra- 
diative feedback there is a clear increase of the median stellar mass 
and a huge decrease in the fraction of brown dwarfs. Thus, in agree- 
ment with earlier works (Bate 2009c, Offner et al. 2009, Urban 
|et al . 2010, Krumholz et al. 20111, we find radiative feedback is 
crucial for determining the IMF even for low-mass star formation. 

Note that, in fact, the calculation produces a protostellar mass 
function (PMF) rather than an IMF ( [Fletcher & Stahler||1994a|b| 
[McKee & Offlier|[20T0) because some of the objects are still ac- 
creting when the calculation is stopped. In this paper, we refer to 
the mass function as an 'IMF' because we compare it to the ob- 
served IMF since the PMF cannot yet be determined observation- 
ally. However, it should be noted that how a PMF transforms into 
the IMF depends on the accretion history of the protostars and how 
the star formation process is terminated. One issue that can be stud- 
ied from the calculation presented here is whether the distribution 
of stellar masses evolves in form during the formation of the stel- 
lar cluster or whether the mass distribution is 'fully-formed' so that 



no matter when the distribution is examined it is always consis- 
tent with being drawn from a constant underlying mass function. 
From Fig. [3] the median stellar mass and the overall shape of the 
distribution obtained using a barotropic equation of state does not 
change with time, except that the maximum mass increases with 
time. This is in contrast to the isothermal calculation of lKrumholzl 
|et al.|(20TT| , who find that the median mass increases by a factor of 
two during their calculation without radiative feedback. This differ- 
ence is probably due to the very different initial conditions, but the 
fact their sink particles are 800 times larger than those used here so 
that they cannot capture small-scale fragmentation may also play a 
role. 

In Fig. |6] we plot the cumulative stellar mass distributions at 
four different times during the radiation hydrodynamical calcula- 
tion. Performing Kolmogorov-Smirnov tests comparing the final 
distribution to each of the early distributions shows that each inter- 
mediate distribution is consistent with being drawn randomly from 
the same distribution as the final distribution. Of course the inter- 
mediate and final distributions are not statistically independent, but 
the test still shows that the stellar mass distribution keeps the same 
form as it evolves during the formation of the cluster This also 
means that in stopping the calculation at t = 1.20 ts we do not 
seem to have stopped at a special point in the evolution, at least as 
far as the mass function is concerned. The only thing that changes 
is that as the number of stars increases with time, so the maximum 
stellar mass increases. Again, this is in contrast to Krumh olz et al.| 
(j20TTjl, who find that with radiative feedback their median stellar 
mass increases by almost a factor of 4 from 0.55 M© when 10% 
of the cloud's mass has been converted to stars to 2 M© when 50% 
the stars contain 50% of the total mass. As mentioned above, this 
difference is probably due to a combination of the denser, centrally- 
condensed initial conditions, the radiative feedback from sink par- 
ticles, and the use of much larger sink particles. 

The differential form of the IMF at the end of the radiation hy- 
drodynamical calculation is shown in the right-hand panel of Fig.|7] 
and is compared with the parameterisations of the observed IMF 
given by |Chabrier| pOOS} , [Kroupal ( |200T] l, and |Salpeter| ( |T955| l. In 
the left-hand panel of the figure, we provide the IMF from the main 
barotropic calculation ofB ate| ( [2009a[ l at the same time for compar- 
ison. In agreement with the smaller radiation hydrodynamical cal- 
culations of |Bate|j2009c| l, the introduction of the radiative feedback 
has clearly addressed the problem of the over-production of brown 
dwarfs and low-mass stars that occurs when using a barotropic 
equation of state (Bate 2009a). In fact, comparing the histogram of 
objects with the parameterisation of the observed IMF by Chabrier| 
( |2005^ , the agreement is almost too good to be true. The cumu- 
lative mass distribution from the calculation (solid line) is com- 
pared with that of |Chabrier| ( |2005) (dot-long-dashed line) in Fig.|5] 
A Kolmogorov-Smimov test gives a 50.5 percent probability that 
the radiation hydrodynamical IMF could have been drawn from 
|Chabrier] s fit to the observed IMF (i.e. the two mass functions are 
indistinguishable) . 

However, despite more than a decade of observational work, 
the form of the IMF in the sub-stellar regime is still quite uncer- 
tain. Although it is now generally accepted that the number of 
stars is larger than the number of brown dwarfs in Galactic star- 
forming regions ( |Chabrier|2003||Greissl et al.|2007[|Andersen et al.| 
[2008) ), considerable uncertainty remains. Rather than trying to de- 
termine the exact form or slope of the substellar IMF, a popular 
method is to compare the number of brown dwarfs in an observed 
region to the number of stars with masses less than that of the 
Sun. [Andersen et ar]j2008^ find that the ratio of stars with masses 
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Figure 7. Histograms giving the initial mass function of the 590 stars and brown dwarfs at t = 1.20tfi from the main barotropic calculation of Bate (2009a) 
(left), and the 183 objects formed at the same time in the radiation hydrodynamical calculation (right). The double hatched histograms are used to denote those 
objects that have stopped accreting, while those objects that are still accreting are plotted using single hatching. The radiation hydrodynamical calculation 
produces far fewer brown dwarfs and low-mass stars and more stars with masses 1.5 Mq and is in good agreement with the Chabiier (2005) fit to the 
observed IMF for individual objects. Two other parameterisations of the IMF are also plotted: Salpeter (1955) and Kroupa (2001). 



0.08 - 1.0 M0 to brown dwarfs with masses 0.03 - 0.08 Mq 
is N{Om - 1.0)/iV(0.03 - 0.08) « 5 ± 2. By combining the 
results of two radiation hydrodynamical calculations of star for- 
mation in 5O-M0 molecular clouds, |Bate| |2009c| found a ratio 
of stars to brown dwarfs of 25/5 « 5. This number is in agree- 
ment with observational constraints, but the statistical uncertainty 
is large. Here we obtain a ratio of 7V(0.08- 1 .0)/Af (0.03-0.08) = 
117/31 ^ 3.8. Eight of theSl low-mass objects and 84 of the 117 
stars were still accreting when the calculation was stopped, so there 
is some uncertainty in this figure due to unknown future evolution. 
But the value is well within observational uncertainties. For the 
main barotropic calculation of Bate ( 2009aj , this ratio was much 
lower: 212/146 « 1.45 at t = 1.20 ts and 408/326 = 1.25 at 
t = 1.50 ts- 

Below 0.03 Mq, the IMF is very poorly constrained, both ob- 
servationally and from the calculation presented here. The radia- 
tion hydrodynamical calculation produced 6 objects with masses 
in this range, with a minimum mass of 17.6 Mj. All of these ob- 
jects had stopped accreting when the calculation was stopped. This 
is very different to the main barotropic calculation. At the same 
time (1.20 te), the barotropic calculation had produced 217 brown 
dwarfs (40 still accreting) with masses less than 30 Mj with a min- 
imum mass of only 5.6 Mj. Even discounting objects that were 
still accreting, the inclusion of radiative feedback has cut the pro- 
duction of these very-low-mass brown dwarfs by a factor of ~ 30 
and significantly increased the minimum mass. It is interesting to 
note that the minimum mass is substantially higher than the opac- 
ity limit for fragmentation 1 Low & Lynden-Bell|1976[|Rees|1976] 
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|Silk|1977a|b||Boyd & Whitworth|2005| l. This is because the opacity 
limit provides a minimum mass, but generally objects will accrete 
from their surroundings to greater masses. Perhaps more impor- 
tantly, the minimum mass is also greater than the estimated masses 
of some objects observed in star-for ming regions ([Z apatero Os orio| 
et al. 2000, Kirkpatrick et al."2001', 'Zapatero Osorio et al. 20W^ 
kirkpatrick et a l. 2006; Lodieu et al. 2008 Luhma n et al. 2008 
2009[ [Weights et al.((2009; Burge ss et al.||2009 
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Luhman et al.||2009| [Quanz et al.||2010| l. While an exact cut off 
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Figure 8. Time of formation and mass of each star and brown dwarf at the 
end of the radiation hydrodynamical calculation. It is clear that the objects 
that are the most massive at the end of the calculation are actually some of 
the first sink particles to form. Objects that are still accreting significantly 
at the end of the calculation are represented with vertical arrows. The hori- 
zontal dashed line marks the star/brown dwarf boundary. Time is measured 
from the beginning of the calculation in terms of the free-fall time of the 
initial cloud (top) or years (bottom). 



presented here do imply that brown dwarfs with masses H, 15 Mj 
should be very rare. 



5.2.7 The origin of the initial mass function 



is difficult to determine from either numerical simulations or ob- 
servations, the results of the radiation hydrodynamical calculation 



[Bate & Bonnell| j2005 ^ analysed two barotropic star cluster for- 
mation simulations that began with clouds of different densities to 
determine the origin of the IMF in those calculations (see also jBate] 
|2005| l. They found that the IMF resulted from competition between 
accretion and 'ejection'. There was no significant dependence of 
the mean accretion rate of an object on its final mass. Rather, there 
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Figure 9. The time-averaged accretion rates of the objects formed in the 
radiation hydrodynamical calculation versus their final masses. The accre- 
tion rates are calculated as the final mass of an object divided by the time 
between its formation and the termination of its accretion or the end of the 
calculation. Objects that are still accreting significantly at the end of the cal- 
culation are represented with horizontal arrows. There is no dependence of 
mean accretion rate on final mass for objects with less than ~ 0.5 Mq (just 
a dispersion). However, there is a low-accretion rate region of exclusion 
for the most massive objects since only objects with mean accretion rates 
greater than their mass divided by their age can reach these high masses dur- 
ing the calculation. The horizontal solid line gives the mean of the accretion 
rates: 1.54 X 10~^ Mq yr~^. The accretion rates are given in M0/tff on 
the left-hand axes and Mq yr~ ^ on the right-hand axes. The vertical dashed 
line marks the star/brown dwarf boundary. 



was a roughly linear correlation between an object's final mass and 
the time between its formation and the termination of its accretion. 
Furthermore, the accretion on to an object was usually terminated 
by a dynamical interaction between the object and another system. 
Note that such an interaction does not necessarily require that the 
object is ejected from the cluster. Many times this is the case, but 
moving an object into a lower density part of the cloud (e.g. out of 
its natal core) or substantially increasing the object's speed with- 
out it becoming unbound can also dramatically reduce its accretion 
rate (c.f., the Bondi-Hoyle accretion formula M oc p/(Cg +i>^)^^^, 
where v is the velocity of the object relative to the gas). Thus, |Batel 
|& Bormell] found that objects formed with very low masses (a few 
Jupiter masses) and accreted to higher masses until their accretion 
was terminated, usually, by a dynamical encounter. This combina- 
tion of competitive accretion and stochastic dynamical interactions 
produced the mass distributions, and [Bate & B onnell (2005) pre- 
sented a simple semi-analytic model which could describe the nu- 
merical results in which the characteristic stellar mass was given 
by the product of the typical accretion rate and the typical time 
between a object forming and having a dynamical interaction that 
terminated its accretion. Bate ( 2009a I found the IMF in their larger 
barotropic calculations also originated in this manner. They found 
the mean accretion rate of a low-mass star did not depend on its 
final mass, but that objects that accreted for longer ended up with 
greater masses and that protostellar accretion was usually termi- 
nated by dynamical interactions. Here we analyse the radiation hy- 
drodynamical calculation using the same methods. 

In Fig. [8] we plot the final mass of an object versus the time 
at which it formed (i.e. the time of insertion of a sink particle). It 



Figure 10. The time between the formation of each object and the termina- 
tion of its accretion or the end of the radiation hydrodynamical calculation 
versus its final mass. Objects that are still accreting significantly at the end 
of the calculation are represented with arrows. As in past barotropic calcu- 
lations, there is a clear linear correlation between the time an object spends 
accreting and its final mass. The solid line gives the curve that the objects 
would lie on if each object accreted at the mean of the time-averaged accre- 
tion rates. The accretion times are given in units of the tg on the left-hand 
axes and years on the right-hand axes. The vertical dashed line marks the 
star/brown dwarf boundary. 

is clear that the most massive stars at the end of the calculation 
were some of the first to begin forming. During the calculation, 
as other lower-mass stars have formed and some have had their 
accretion terminated, these stars have continued to grow to higher 
and higher masses. [Maschberger et al.| ( |2010^ have argued that such 
a cluster formation process naturally produces a relation between 
cluster mass and maximum stellar mass similar to that which is 
observed ( [Weidner & Kroupa|2006[[Weidner et al.|2010[ >, although 
others argue that the observations are also consistent with random 
sampling from a universal IMF ( [Lamb et al.|2010] [Fumagalli, da| 
|Silva&Krum holz|2011>. 

In Figs.[9][T0] and|l I[ we plot similar figures to those found in 
[Bate & Bonnell| ( |2005| and |Batelp005l|2009a|b| l. These figures dis- 
play the same trends as found in the barotopic calculations. Fig.|9] 
gives the time-averaged accretion rates of all the objects formed in 
the radiation hydrodynamical calculation versus the object's final 
mass. The time-averaged accretion rate is the object's final mass 
divided by the time between its formation (i.e. the insertion of a 
sink particle) and the end of its accretion (defined as the last time 
its accretion rate drops below lO"'^ M© yr~^) or the end of the 
calculation. As in the barotropic calculations, there is no depen- 
dence of the time-averaged accretion rate on an object's final mass, 
except that objects need to accrete at a rate at least as quickly as 
their final mass divided by their age (i.e., the lower right potion 
of Fig.[9]cannot have any objects lying in it). This means that the 
most massive stars have higher time-averaged accretion rates than 
the bulk of the stars and VLM objects. But, on the other hand, if the 
calculation were continued longer, objects that accrete with lower 
time-averaged accretion rates could also reach high masses. Note 
that these results should not be used to infer that the typical ac- 
cretion rate remains independent of mass above 3 Mq. The calcu- 
lation presented here does not produce any high-mass protostars, 
but other studies have reported that the accretion rates of protostars 
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Figure 11. For each single object that has stopped accreting by the end of 
the main calculation, we plot the time of the ejection of the object from 
a multiple system versus the time at which its accretion is terminated. As 
in past barotropic calculations, these times are correlated, showing that the 
termination of accretion on to an object is usually associated with dynam- 
ical ejection of the object. Open circles give those objects where multiple 
'ejections' are detected by the ejection detection algorithm and, hence, the 
ejection time is ambiguous (see the main text). Binaries have been excluded 
in the plot because it is difficult to determine when a binary has been ejected. 



with masses <; 3 Mq do increase with mass (e.g., [Urban et al.|2010| 
[Krumholz et al.|201 

The mean of the accretion rates is 1.5 x 10~^ yr~^, 
which is within a factor of ~ 2 of the mean accretion rates of 
the barotropic calculations in all of the above papers. Thus, the 
mean accretion rate does not depend significantly on cloud den- 
sity (Bate & Bonnell 2005 ), the equation of state of high-density 
gas (Bate 2005), the total mass of the gas cloud (Bate 2009al, or 
on whether the calculation is performed using a barotropic equa- 
tion of state or radiation hydrodynamics. It only depends on the 
mean temperature of the initial cloud ( [Bate|2005| l, in that it scales 
roughly as T^^^ (or, equivalently, cf/G, where is the mean sound 
speed on large-scales ; |Shu| 1971\ . The dispersion in the accretion 
rates is about 0.37 dex, also similar to the previous barotropic sim- 
ulations. Thus, rather than the final mass of a star depending on its 
average accretion rate, the primary determinant of the final mass of 
a star or brown dwarf is the period over which it accretes. Fig.|10| 
very clearly shows the linear relation (with some dispersion) be- 
tween the period of time over which an object accretes and its final 
mass. This means that the higher characteristic stellar mass pro- 
duced when radiative feedback is included is due to an increase in 
the average time over which an object accretes. 

In Fig. [TT] for each object that has stopped accreting by the 
end of the main calculation (excluding the components of bina- 
ries), we plot the time at which the object undergoes a dynamical 
ejection versus the time that its accretion is terminated. The strong 
correlation shows that accretion is usually terminated by a dynam- 
ical encounter with other objects, as seen in the barotropic calcu- 
lations. We define the time of ejection of an object as the last time 
the magnitude of its acceleration drops below 2000 km s^^ Myr^^ 
( |Bate|[2009a^ or the end of the calculation. The acceleration cri- 
terion is based on the fact that once an object is ejected from a 
stellar multiple system, sub-cluster, or cluster through a dynami- 
cal encounter, its acceleration will drop to a low value. We exclude 
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Figure 12. The initial mass functions produced by the radiation hydrody- 
namical calculation (histogram) and a comparison with the fit using the sim- 
ple accretion/ejection IMF model (thick curve) of |Bate & Bonnell] j2005^ . 
Statistically, the hydrodynamical and the model IMFs are indistinguishable 
(a Kolmogorov-Smirnov test gives a 19 percent probability that the hydro- 
dynamical IMF could have been drawn from the model IMF). Also shown 
are the Salp eter| jl955^ slope (solid straight line), and the |Kroupa| pOOH 
(solid broken line) and [Chabrier (2005) (thin curve) mass functions. The 
vertical dashed line is the stellar-substellar boundary. 



binaries because they have large accelerations throughout the cal- 
culation which frequently results in false detections of ejections. 
Also, in Fig. [TT] we use two different symbols (filled circles and 
open circles). For the former we are confident of the ejection time. 
However, for those objects denoted by the open circles, we find 
that at least two 'ejections' more than 2000 years apart have oc- 
curred. These are usually objects that have had a close dynamical 
encounter with a multiple system that has put them into long-period 
orbits rather than ejecting them. In these cases, we chose the 'ejec- 
tion' time closest to the accretion termination time but we use an 
open symbol to denote our uncertainty in whether or not we have 
identified the best time for the dynamical encounter 

We find that, excluding binaries, for 40 objects out of 47 (85%) 
the accretion termination time and the ejection time are within 2000 
years of each other If we also exclude those objects for which 
we are uncertain in our identifications of the ejection times as de- 
scribed above, we find 33 objects out of 40 (83%) are consistent 
with ejection terminating their accretion. These are probably lower 
limits in the sense that it is difficult to determine in an automated 
way the time at which an ejection occurs and an erroneous value 
is much more likely to differ from the accretion termination time 
by more than 2000 years than coincide with it. In any case, it is 
clear that for the majority of objects their accretion is terminated 
by dynamical encounters with other stellar systems. 

In Fig. [12] we compare the IMF obtained from the radi- 
ation hydrodynamical calculation with the semi-analytic accre- 
tion/ejection IMF model of [Bate & Bonnell] ( |2005^ using parame- 
ters determined from the radiation hydrodynamical calculation (see 
also Bate 2009b l. The parameters are: the mean accretion rate and 
its dispersion (given above), period of time over which stars form 
(i.e. 90,000 yrs), the characteristic ejection time, and the minimum 
stellar mass. The characteristic ejection time, reject = 62, 400 yr, 
is chosen such that the mean number of objects that have finished 
accreting over the time period equals that from the radiation hy- 
drodynamical calculation (64 objects). The minimum stellar mass 
primarily determines the minimum mass cut-off to the IMF, rather 
than the shape of the rest of the IMF. For the semi-analytic IMF in 
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Figure 13. The magnitudes of the velocities of each star and brown dwarf 
relative to the centre-of-mass velocity of the stellar system at the end of the 
radiation hydrodynamical calculation. For binaries, the centre-of-mass ve- 
locity of the binary is given, and the two stars are connected by dotted lines 
and plotted as squares rather than circles. Objects still accreting at the end 
of the calculation are denoted by horizontal arrows. The root mean square 
velocity dispersion for the association (counting each binary once) is 5.5 
km s^i (3-D) or 3.2 km s-^ (1-D). There is a dependence of the velocity 
dispersion on mass with VLM objects having a lower velocity dispersion 
than stars (see the main text). Binaries are found to have a slightly lower 
velocity dispersion than single objects of only 4.6 km s~^ (3-D). The ver- 
tical dashed line marks the star/brown dwarf boundary. 



Fig.[T2]we choose 5 Jupiter-masses, but 10-15 Jupiter-masses result 
in similarly good fits. A Kolmogorov-Smimov test comparing the 
semi-analytic IMF to the IMF obtained from the radiation hydro- 
dynamical calculation shows that the latter is consistent with being 
randomly drawn from the former (probability 19%). 

In conclusion, the origin of the IMF in the radiation hydrody- 
namical calculation is the same as in the past barotropic calcula- 
tions: the IMF originates from competition between accretion and 
dynamical encounters. Objects end up with low masses if their ac- 
cretion is terminated (by a dynamical encounter) soon after they 
form. Objects end up with high masses by accreting for an ex- 
tended period. The reason the characteristic stellar mass is larger 
when radiative feedback is included is that objects typically accrete 
for longer before their accretion is terminated. This is because the 
radiative feedback increases the typical distance between objects 
( |Bate|2009c| >, and so dynamical interactions take longer to occur. 



In Fig. [T3] we plot the magnitude of the velocity of every star 
or brown dwarf relative to the centre of mass of the stellar system 
at the end of the radiation hydrodynamical calculation. For bina- 
ries (including those that are sub-components of triples and quadru- 
ples), we plot the two components with the centre of mass velocity 
of the binary using filled squares connected by a dotted line. The 
overall root mean square (rms) velocity dispersion (counting each 
binary only once) is 5.5 km s~^ (3-D) or 3.2 km s~^ (1-D). This 
is almost identical to the velocity dispersion found by |Bate| l [2009a[ > 
at the end of the main barotropic calculation, and larger than the 
velocity dispersions found from calculations of star formation in 
small 5O-M0 clouds ( |Bate et al. 2003 ; Bate & Boniiei l 2005 ; Bate] 
|2005, 2009b I. The velocity dispersion is 40% higher than the ID 
velocity dispersion of 2.3 km s~^ seen in the Orion Nebula Cluster 
l |Jones & Walkerl 19881 |Tian et al.|I996| l, consistent with the fact 
that the system formed here is lower in mass, but denser. 

[Reipurth & Clarke] ( [200 1 [ l suggested that a greater velocity 
dispersion for brown dwarfs than stars may be a possible signa- 
ture that brown dwarfs form as ejected stellar embryos. Past A'^- 
body simulations of the breakup of small-A'^ clusters with A'^ > 3 
l |Sterzik & Durise n"I998") and SPH calculations of TV = 5 clus- 
ters embedded in gas ( Delgado-Donate et al. 2003 1 found that there 
was no strong dependence of the velocity of an object on its mass, 
but both found that binaries should have a smaller velocity disper- 
sion than single objects due to the recoil velocities of binaries being 
lower. On the other hand, Delga do-Donate et al.| ( |2004 ) performed 
simulations of star formation in small turbulent clouds and found 
that the velocity dispersions of singles and binaries were indistin- 
guishable, but that higher-order multiples had significantly lower 
velocity dispersions. 

From the large barotropic calculations of 'Bate (2009a), it was 
found that stars tend to have a slightly higher dispersion than VLM 
objects and that binaries have a lower velocity dispersion than sin- 
gle objects. These same relations are also found from the radiation 
hydrodynamical calculation. The rms velocity dispersion of VLM 
systems is 4.1kms^^ (3-D) while for the stars (masses^ O.IMq) 
the rms velocity dispersion is 6.5 km s^^ (3-D). Binaries (most 
of which have stellar primaries) have a velocity dispersion of 4.6 
km s~^ (3-D), lower than both the velocity dispersion of all stars 
and the overall velocity dispersion. 

Observationally, while there is no strong evidence for VLM 
objects having a different velocity dispersion to stars, or binaries 
having a different velocity dispersion to single objects, studies of 
the radial velocities of stars and brown dwarfs in the Chamaeleon I 
dark cloud do find that brown dwarfs have a marginally lower ve- 
locity dispersion than the T Tauri stars ( |Joergens" & Guenther|200Tj 
|Joergens|2006[ l, in qualitative agreement with the results from both 
the barotropic calculations of |Bate|f2009a| l and radiation hydrody- 
namical calculation discussed here. 



3.3 Stellar kinematics 

At the end of the main barotropic calculation, most of the sub- 
clusters had merged together and the stellar distribution essentially 
consisted of a single cluster surrounded by a halo of ejected stars 
and brown dwarfs. |Bate| j2009a| l analysed the radial properties of 
this cluster, looking for evidence of mass segregation and for radial 
variations of the stellar velocity dispersion and binarity. Unfortu- 
nately, we are not able to follow the radiation hydrodynamical cal- 
culation until all the sub-clusters and filaments have collapsed into 
a single cluster and we have a lot fewer stars and brown dwarfs, so 
we restrict ourselves only to discussing the kinematics of the pop- 
ulation as a whole without attempting to look for spatial variations. 



3.4 Stellar encounters and disc sizes 



Reip urth & Clarke] p001| l also speculated that if brown dwarfs 
formed via ejection, they might have smaller, lower-mass discs 
than stars. As mentioned in Section [3.2.1| the IMF in the radiation 
hydrodynamical calculation presented here and in the many past 
barotropic star cluster formation calculations originates through 
competition between accretion and ejection, but this applies both to 
stars and brown dwarfs. The only difference is that brown dwarfs 
are ejected soon after they form, before they have accreted much 
mass, while stellar ejections occur after a longer period of accre- 
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Figure 14. The closest encounter distances of eacli star or brown dwarf dur- 
ing tlie radiation iiydrodynamical calculation versus the final mass of each 
object. Objects that are still accreting significantly at the end of the calcu- 
lation are denoted with arrows indicating that they are still evolving and 
that their masses are lower limits. Binaries are plotted with the two compo- 
nents connected by dotted lines and squares are used as opposed to circles. 
The horizontal dotted line marks the size of the sink particle accretion radii, 
within which there is no gas. The vertical dashed line marks the star/brown 
dwarf boundary. Note that two objects almost merged with one another dur- 
ing the calculation (the merger radius was set to be 0.01 AU). 



Figure 15. Due to dynamical interactions, stars and brown dwarfs poten- 
tially have their discs truncated to approximately 1/2 of the periastron sepa- 
ration during the encounter (see also Figure [T4) . At the end of the radiation 
hydrodynamical calculation, we plot the cumulative fraction objects as a 
function of the potential truncation radius. We exclude binaries and any ob- 
jects that are still accreting at the end of the calculation. The solid line gives 
the result for all stai's and brown dwarfs, while the dotted, short-dashed, 
and long-dashed lines give the cumulative distributions for the mass ranges 
M < 0.1, 0.1 s£ M < 0.3, and 0.3 ^ M < 1.0, respectively. There 
are no stars with masses M ^ 1.0 Mq that are not either accreting or in 
multiple systems at the end of the calculation. 



tion. Discs around both stars and brown dwarfs may be truncated 
by dynamical encounters and ejections. 

In Figure[T4] we plot the distance of the closest encounter that 
every star and brown dwarf has had by the end of the radiation hy- 
drodynamical calculation. As in past barotropic calculations, there 
is a wide range of closest encounter distances (including two ob- 
jects that almost merged), but the closest encounters tend to have 
occurred for stars rather than brown dwarfs. Dynamical encounters 
between objects will truncate any circumstellar discs. However, this 
plot cannot be taken to mean that many stars have small discs be- 
cause of several reasons. First, as will be seen in Section [331 mul- 
tiplicity is a strong function of primary mass. In Figure[T4]it clear 
that binaries are responsible for many of the 'closest encounters'. 
Second, objects that are still accreting at the end of the calculation 
are still evolving and, since the mass of an object depends on its 
'age' more massive accreting objects are more likely to have had 
close encounters. Finally, as noted by |Bate et al.H2003] >, many stars 
that have close encounters have new discs formed from accretion 
subsequent to their closest dynamical encounter. 

Despite these difficulties, if an object suffers a dynamical en- 
counter that terminates its accretion this encounter will truncate any 
disc that is larger than approximately 1/2 of the periastron distance 
during the encounter ( |HaIl, Clarke & Pringle|1996) . Therefore, ex- 
cluding binaries and objects that are still accreting, determining the 
distribution of 1/2 of the closest encounter distance should give us 
an indication of the disc size distribution around single objects that 
have reached their final masses. Note that formally we have still 
included the wide components of triple and quadruple systems, but 
these constitute only 13 objects out of the 94 'single' non-accreting 
objects and all but one are in orbits with semi-major axes greater 
than 10 AU so their inclusion should not adversely affect our con- 
clusions. 



In Figure[T5] we plot the cumulative distributions of disc trun- 
cation radii (taken to be 1/2 of the closest encounter distance) for 
these objects. The solid line gives the cumulative distribution for 
all 94 objects, while in the other distributions we break the sample 
into mass bins of Af, < 0.1 M©, 0.1 M© ^ M, < 0.3 M©, 
and 0.3 Mq sC M. < 1.0 Mq . [Bate| p009al l found that in the 
main barotropic calculation there was a clear relation such that 
more massive stars tended to have had closer encounters. Here, the 
statistics are not as good, but the stars in the highest mass bin have 
clearly typically had closer encounters. 

For VLM objects, 20% have truncation radii greater than 40 
AU, while 1/2 have truncation radii greater than 10 AU. It has been 
known for a decade from infrared excess that young brown dwarfs 
have discs (|Muench et al.||2001| |Natta & Testi||200T] |Apai et al 



2002) |Natta et aI.|2002[|Jayawardhana et al.|2003[|Luhman et al. 



|2005 |Moni n et al.'20I0), some of which also display evidence for 
accretion jjayawardhana et al. 2002). At least some of these discs 
are inferred to have radii of 20-40 AU (e.g. [Luhman et al.|2007] l, 
while [Scholz et al^p006[ l estimate that at least 25% of the brown 
dwarfs they survey in Taurus have discs larger than 10 AU in ra- 
dius. The cumulative distribution of truncation radii in Fig. |15|is 
consistent with these observations, but it should be used with cau- 
tion. First, the simulation presented here produces a dense stellar 
cluster Disc truncation may be less important for setting disc sizes 
in a lower-density star-forming region like Taurus. Second, Figure 
[15] does not give a disc size distribution. At best, it is a distribu- 
tion of lower limits to disc sizes because of the fact that stars can 
suffer a close dynamical encounter, but then accrete more mate- 
rial from the molecular cloud and form a new disc. This happens 
frequently in the simulation, especially for the higher mass stars. 
However, the distribution may be more useful for VLM objects be- 
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Figure 16. We find the closest object to each star or brown dwarf when it 
forms and plot their final versus initial separation (open circles). We also 
plot the final semi-major axes versus the initial separations of all binaries at 
the end of the calculation (small filled circles). Note that the closest object 
when a star or brown dwarf forms often does not remain close. Also, many 
of the close multiple systems at the end of the calculation are composed of 
objects that formed at large distances from each other. These results indi- 
cate the importance of dynamical interactions and orbital decay during the 
calculation. 



cause they tend to have their accretion terminated soon after they 
form by dynamical encounters and generally will not subsequently 
accrete significantly from the molecular cloud. 



there exist 2 1 binary systems and one triple system with separations 
< 10 AU. 

The mechanisms by which close binaries form have been dis- 
cussed in detail by |Bate et al. ( 2002b) . Briefly, Bate et al. found that, 
rather than forming directly via fragmentation, the close binary sys- 
tems form from the orbital decay of wider systems through a com- 
bination of dynamical interactions, accretion, and the interaction of 
binaries and triples with circumbinary and circumtriple discs. Dy- 
namical interactions can harden existing wide binaries by removing 
angular momentum and energy from their orbits. They also pro- 
duce exchange interactions in which a temporary unstable multiple 
system decays by ejecting one of the components of the original 
binary (usually the lowest-mass object). However, dynamical in- 
teractions alone cannot produce the observed frequency of close 
binary systems, either beginning with stellar clusters ([Kroupa & 



Burkert'200r) or during the dissolution of small-N clusters |Sterzik 



& Durisen 19981. The key is to have dissipative dynamical inter- 
actions, where the presence of gas allows dynamical encounters 
to dissipate energy and transport angular momentu m ^Bate et al.| 
2002b}. A good example is a star-disc encounter dLarson 1990 



parke&Pringle'1991a''b', ' McDonald & Clarke|1995||Heller|1995 
Hall, Clarke & Pringle 1996| . In addition to dynamical interactions. 



accretion onto a binary from an envelope decreases the binary's 
separation unless the specific angular momentum of the accreted 
material is significantly greater than that of the bina ry (| Arty mow-] 
[iczl[T983l |Batel[T997l [Bate & BonneIil[T997l [Bate 2000). It can 
also change the relative separations of a triple system, destabilis- 
ing it and forcing dynamical interactions ( [Smith, Bonnell & Bate| 
|1997[ ). Circumbinary discs can remove large amounts of orbital 
angular momentum from an embedded binary system via gravita- 
tional torques, thus tightening its orbit jPringle| 1 99 1 1 1 Artymowicz] 
|etal.|1991i l. 

These mechanisms work in the calculation here to produce 
the close systems. They also reduce the separations of wide ~ 
1000 AU systems to systems with intermediate separations ~ 10 — 
100 AU. 



3.5 The formation of multiple systems 

The opacity limit for fragmentation sets a minimum initial binary 
separation of ~ 10 AU since the size of a slowly-rotating first 
hydrostatic core is ~ 5 AU fLarso n|1969[ [Masunaga & Inutsuka| 
[2500||Tomida et al.|2010 : Bate 20li|l. In Fig.|16| we plot the dis- 
tance to the closest other star or brown dwarf when each star or 
brown dwarf forms and the distance to this object at the end of the 
calculation (open circles). Also plotted is the initial separation and 
final semi-major axis of all binaries at the end of the calculation 
(filled circles). No object forms closer than 10 AU from an exist- 
ing object, consistent with the expectations from opacity limited for 
fragmentation. 

Examining Fig.[T6] we find that some objects that begin with 
close separations end up well separated. Such situations occur 
when one of the objects is involved in a dynamical interaction (e.g. 
ejection from a group or multiple system). Alternately, objects that 
are initially widely separated (100 — 10* AU) can end up in close 
bound systems (separations < 100 AU). In fact, most of the close 
binaries at the end of the calculation are composed of mutual near- 
est neighbours at the time of formation (i.e. in the figure the filled 
circles are inside the open circles), but the separations of these ob- 
jects may have been reduced by up to 3 orders of magnitude during 
the evolution. In particular, despite the fact that no objects form 
closer than « 10 AU from each other, at the end of the calculation 



3.6 Multiplicity as a function of primary mass 

We turn now to the properties of the binary and higher-order mul- 
tiple stars and brown dwarfs produced by the simulation. Observa- 
tionally, it is clear that the fraction of stars or brown dwarfs that 
are in multiple systems increases with stellar mass (massive stars: 



Mason et allfT998 Jreibisch et al.|[T999i [Shatsky & Tokovinin 



2002^ 



Kobulnick y & Fryer|2007[ [Kouwenhoven et^al.|2007 



Mason 



et al.|2009i intermediate-ma ss stars: 'Patience et al.'2002', solar-type 



stars :|Duquennoy & Mayor, 1 991 , Raghavan et al. 2010, M-dwarfs: 
Fischer & Marcy 1992 and very-low-mass stars and brown dwarfs: 



Close et al.,.2003, .Siegler et al.|[2005] [Basri & Reiriersl[2006l >. It 
also seems that the multiplicity of young stars in low-density star- 
forming regions is somewhat higher than that of field stars ('Leinert 
|et al. 1993, Ghez et al. 1993; Simon et al. 1995; Du chene et al 



2007 ). However, IC348 has a similar binary frequency to the field 
(Du chene et al.|1999| >. In the Orion Nebula Cluster, [Kbhler et al.| 
(2006} find that the binary frequency of low-mass stars is similar to 
that of field M dwarfs and lower than that of field solar-type stars, 
but that stars with masses A/ > 2 M0 have a higher binarity than 
stars with 0.1 < M < 2 M0 by a factor of 2.4 to 4. 

To quantify the fraction of stars and brown dwarfs that are in 
multiple systems we use the multiplicity fraction, m/, defined as a 
function of stellar mass. We define this as 
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Object Numbers Mi M2 q a e Initial Relative Spin Spini -Orbit Spin2-Orbit Comments 

Separation or Orbit Angle Angle Angle 
[M0] [M0] [AU] [AU] [deg] [deg] [deg] 



138, 117 


0.29 


0.29 


0.97 


0.39 


0.61 


907 


91 


129 


64 




32, 50 


1.86 


0.59 


0.32 


0.88 


0.73 


748 


105 


68 


-94 




20, 23 


2.76 


1.97 


0.72 


0.88 


0.76 


1533 


13 


49 


62 




25, 26 


1.81 


1.23 


0.68 


0.91 


0.61 


419 


7 


31 


27 




38, 45 


1.40 


1.20 


0.86 


1.08 


0.60 


88.3 


4 


69 


73 




6, 13 


1.90 


1.30 


0.68 


1.20 


0.71 


950 


28 


55 


37 




1 7 


1.37 


1.07 


0.78 


1.25 


0.63 


366 


87 


96 


-169 




64, 79 


0.84 


0.10 


0.12 


1.95 


0.24 


100 


48 


28 


57 




35 33 


1.91 


1.50 


0.79 


1.96 


0.67 


2833 


88 


97 


77 




87, 121 


1.26 


0.40 


0.32 


2.02 


0.26 


2159 


73 


42 


-36 




118 163 


0.59 


0.21 


0.35 


2.24 


0.16 


210 


25 


43 


66 




3, 8 


2.27 


1.00 


0.44 


2.56 


0.45 


3268 


21 


13 


17 




98 109 


0.29 


0.21 


0.71 


3.01 


0.44 


25.4 


4 


68 


66 




5, 16 


2.53 


2.49 


0.98 


3.27 


0.26 


57.3 


3 


23 


23 




116, 132 


0.25 


0.12 


0.45 


4.82 


0.52 


547 


19 


59 


56 




63, 154 


0.29 


0.14 


0.50 


6.44 


0.62 


22.1 


6 


1 


.5 




115 183 


0.52 


0.04 


0.08 


8.40 


0.45 


209 


32 


12 


42 


*\tarA/T M Irinle 


90, 103 


0.18 


0.10 


0.54 


8.75 


0.10 


63.8 


10 


12 


10 


vStarA^l M binarv 


122, 145 


0.23 


0.14 


0.62 


8.81 


0.23 


158 


13 


21 


13 




102 110 


0.27 


0.26 


0.96 


9.36 


0.12 


67.8 


7 


37 


29 




159 150 


0.16 


0.15 


0.96 


9.91 


0.35 


1138 


13 


38 


51 




59, 68 


0.08 


0.05 


0.61 


10.6 


0.46 


12.6 


13 


17 


4 


VLM binary 


19 27 


0.28 


0.25 


0.89 


12.8 


0.08 


58.4 


3 


22 


22 




76, 83 


0.25 


0.21 


0.84 


12.8 


0.42 


267 


9 


30 


22 




164 179 


0.15 


0.03 


0.21 


14.3 


0.80 


23.3 


22 


4 


22 


Star A/7 M trinle 


92, 133 


0.27 


0.20 


0.74 


14.3 


0.16 


52.6 


20 


19 


-2 




44 82 


1.03 


0.91 


0.88 


14.3 


0.01 


267 


g 


35 


31 




41, 89 


0.25 


0.18 


0.71 


14.7 


0.09 


324 


11 


31 


20 




94, 129 


0.18 


0.09 


0.49 


14.8 


0.09 


445 


56 


57 


3 


StarAa.M binary 


4, 84 


1.33 


1.06 


0.80 


19.3 


0.02 


117 


12 


40 


41 




160, 168 


0.19 


0.15 


0.81 


19.5 


0.45 


902 


15 


29 


34 




104, 93 


0.35 


0.34 


1.00 


22.5 


0.13 


193 


6 


36 


39 




140, 147 


0.09 


0.09 


0.94 


26.1 


0.01 


48.7 


5 


23 


19 


VLM binary 


65, 77 


0.27 


0.26 


0.95 


27.5 


0.01 


396 


7 


21 


14 




88, 127 


0.21 


0.10 


0.50 


29.7 


0.05 


443 


41 


24 


-20 




175, 181 


0.09 


0.08 


0.98 


36.4 


0.06 


64.8 


106 


53 


-54 


VLM binary 


10, 17 


3.84 


0.70 


0.18 


165 


0.16 


6244 


67 


104 


92 




67, 74 


0.14 


0.05 


0.38 


356 


0.35 


327 


74 


144 


77 


StarA'LM binary 


49, 96 


0.38 


0.27 


0.71 


406 


0.74 


969 


44 


70 


46 




106, 148 


0.71 


0.12 


0.17 


474 


0.70 


369 


53 


70 


26 




12, 105 


0.88 


0.04 


0.04 


620 


0.42 


1185 


77 


91 


93 


Star/VLM binary 


153, 182 


0.13 


0.03 


0.26 


721 


0.82 


1289 


44 


128 


93 


Star/VLM binary 


171, 174 


0.10 


0.08 


0.79 


8366 


0.90 


1743 


94 


36 


99 


StarA'LM binary 


(25, 26), 37 


f3 04) 


1.68 


0.55 


5.53 


0.19 




34 








(64, 79), 55 


(0 94) 


0.86 


0.91 


18.1 


0.10 




4 








(115, 183), 149 


(0.56) 


0.19 


0.34 


19.9 


0.10 




14 






StarA'LM triple 


(5, 16), 15 


(5.02) 


2.95 


0.59 


23.5 


0.16 




36 








( 87, 121), 100 


(1.66) 


1.21 


0.73 


36.8 


0.24 




32 








(122, 145), 123 


(0.37) 


0.26 


0.70 


45.2 


0.08 




5 








(116, 132), 131 


(0.37) 


0.15 


0.41 


57.8 


0.19 




6 








(104, 93), 134 


(0.69) 


0.22 


0.31 


108 


0.17 




10 








(164, 179), 173 


(0.18) 


0.07 


0.40 


194 


0.49 




126 






Star/VLM triple 


((87, 121), 100), 139 


(2.87) 


0.06 


0.02 


138 


0.27 












(4, 84), (44, 82) 


(2.39) 


(1.94) 


0.81 


139 


0.03 












(38, 45), (32, 50) 


(2.59) 


(2.45) 


0.94 


142 


0.39 












((115, 183), 149), 126 


(0.75) 


0.50 


0.67 


145 


0.09 












(76, 83), (41, 89) 


(0.46) 


(0.43) 


0.93 


161 


0.45 












((25, 26), 37), 40 


(4.72) 


3.38 


0.72 


177 


0.31 












((116, 132), 131), 119 


(0.52) 


0.11 


0.21 


10575 


0.92 













Table 3. The properties of the 40 multiple systems at the end of the calculation. The structure of each system is described using a binary hierarchy. For each 
'binary' we give the masses of the primary Mi and secondary M2, the mass ratio q = M2/M1, the semi-major axis a, the eccentricity e, the separation of 
the components when the second one first formed, the relative spin angle, and the angles between orbit and each of the primary's and secondary's spins. For 
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Figure 17. Multiplicity fraction as a function of primary mass. The left panel gives the result at the end of the radiation hydrodynamical calculation. On the 
right, we give the result from the main barotropic calculation of Bate (2009a) at the same time. The blue filled squares surrounded by shaded regions give the 
results from the calculations with their statistical uncertainties. The thick solid lines give the continuous multiplicity fractions from the calculations computed 
using a boxcar average. The open black squares with error bars and/or upper/lower limits give the observed multiplicity fractions from the surveys of Close 
et al. (2003), Basri & Reiners (2006), Fisher & Marcy (1992), Raghavan et al. (2010), Duquennoy & Mayor (1991), Preibisch et al. (1999) and Mason et al. 
(1998), from left to right. Note that the error bars of the Duquennoy & Mayor (1991) results have been plotted using dashed lines since this survey has been 
superseded by Raghavan et al. (2010). The observed trend of increasing multiplicity with primary mass is well reproduced by both calculations. Note that 
because the multiplicity is a steep function of primary mass it is important to ensure that similar mass ranges are used when comparing the simulation with 
observations. 



Mass Range [Mq] 


Single 


Binary 


Triple 


Quadruple 


M < 0.03 


7 











0.03 ^ M < 0.07 


20 











0.07 ^ M < 0.10 


8 


3 








0.10 ^ M < 0.20 


17 


7 


1 





0.20 ^ M < 0.50 


21 


9 


2 


2 


0.50 ^ M < 0.80 


5 


2 





1 


0.80 ^ M < 1.2 


2 


1 


1 





M > 1.2 


4 


6 


1 


4 


All masses 


84 


28 


5 


7 



Table 2. The numbers of single and multiple systems for different primary 
mass ranges at the end of the radiation hydrodynamical calculation. 



mf = 



B + T + i 



S + B + T + Q' 



(2) 



where S is the number of single stars within a given mass range 
and, B, T, and Q are the numbers of binary, triple, and quadruple 
systems, respectively, for which the primary has a mass in the same 
mass range. Note that this differs from the companion star fraction, 
csf, that is also often used and where the numerator has the form 
B + 2T + 3Q. We choose the multiplicity fraction following |Hub-| 
|ber & Whitworth ( 2005 ) who point out that this measure is more 
robust observationally in the sense that if a new member of a multi- 
ple system is found (e.g. a binary is found to be a triple) the quantity 
remains unchanged. We also note that it is more robust for simula- 
tions too in the sense that if a high-order system decays because it 
is unstable the numerator only changes if a quadruple decays into 
two binaries (which is quite rare). Furthermore, if the denominator 
is much larger than the numerator (e.g. for brown dwarfs where the 
multiplicity fraction is low) the production of a few single objects 
does not result in a large change to the value of mf. This is useful 
because many of the systems in existence at the end of the calcula- 
tions presented here may undergo further dynamical evolution. By 
using the multiplicity fraction our statistics are less sensitive to this 
later evolution. 



The method we use for identifying multiple systems is the 
same as that used by |Bate| j2009a| l, and a full description of the 
algorithm is given in the method section of that paper. When 
analysing the simulations, some subtleties arise. For example, 
many 'binaries' are in fact members of triple or quadruple systems 
and some 'triple' systems are components of quadruple or higher- 
order systems. From this point on, unless otherwise stated, we de- 
fine the numbers of multiple systems as follows. The number of bi- 
naries excludes those that are components of triples or quadruples. 
The number of triples excludes those that are members of quadru- 
ples. However, higher order systems are ignored (e.g. a quintuple 
system may consist of a triple and a binary in orbit around each 
other, but this would be counted as one binary and one triple). We 
choose quadruple systems as a convenient point to stop as it is likely 
that most higher order systems will not be stable in the long-term 
and would decay if the cluster was evolved for many millions of 
years. The numbers of single and multiple stars produced by the 
radiation hydrodynamical calculation are given in Table [2] follow- 
ing these definitions. In Table [3] we give the properties of the 40 
multiple systems. 

In the left panel of Fig.[T7] we compare the multiplicity frac- 
tion of the stars and brown dwarfs as a function of stellar mass 
obtained from the radiation hydrodynamical calculation with ob- 
servations. The results from a variety of observational surveys (see 
the figure caption) are plotted using black open boxes with associ- 
ated error bars and/or upper/lower limits. The data point from the 
survey of |Duquennoy & Mayor| ( [T99I| l is plotted using dashed lines 
for the error bars since this survey has been recently superseded by 
that of Ra ghavan et al.| ( (20I0| l. The results from the radiation hy- 
drodynamical simulation have been plotted in two ways. First, us- 
ing the numbers given in Table|2]we compute the multiplicity in six 
mass ranges: low-mass brown dwarfs (masses < 0.03 Mq), VLM 
objects excluding the low-mass brown dwarfs (masses 0.03 — 0.10 
Mq), low-mass M-dwarfs (masses 0.10 — 0.20 Mq), high-mass 
M-dwarfs (masses 0.20 — 0.50 Mq), K-dwarfs and solar-type stars 
(masses 0.50 — 1.20 Mq), and intermediate mass stars (masses 
> 1.2 Mq). The filled blue squares give the multiplicity fractions 
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in these mass ranges, while the surrounding blue hatched regions 
give the range in stellar masses over which the fraction is calcu- 
lated and the Icr (68%) uncertainty on the multiplicity fraction. In 
addition, a thick solid line gives the continuous multiplicity fraction 
computed using a boxcar average of the results from the radiation 
hydrodynamical simulation. The width of the boxcar average is one 
order of magnitude in stellar mass. 

The radiation hydrodynamical calculation clearly produces a 
multiplicity fraction that strongly increases with increasing primary 
mass. Furthermore, the values in each mass range are in agree- 
ment with observation. In the right panel of Fig. [TT] we provide 
the equivalent quantities obtained from the main barotropic calcu- 
lation of|^Bate (2009a) at the same time as the end of the radiation 
hydrodynamical calculation. Those readers who wish to examine 
the multiplicity at the end of the barotropic calculation can find this 
in |Bate| ( |2009a^ . The barotropic calculations also produce a mul- 
tiplicity that is a strongly increasing function of mass. In fact, the 
results using radiation hydrodynamics and a barotopic equation of 
state are very similar The main barotropic calculation gives mul- 
tiplicities that are somewhat higher for primary masses > 0.2 
than those given by the radiation hydrodynamical calculation, but 
the results are consistent within the statistical uncertainties. 

It is important to note that the surveys with which we are 
comparing the multiplicities are primarily of field stars rather than 
young stars. This is necessary because surveys of young stars ei- 
ther do not sample a large range of separations and mass ratios, 
or the statistics are too poor However, there may be considerable 
evolution of the multiplicities between the age of the stars when 
the calculations were stopped (~ 10^ yrs) and a field population. 
This question of the subsequent evolution of the clusters produced 
by hydrodynamical simulations was recently tackled by |Moeckel &| 
|Bate| j2010) who took the end point of the main barotropic calcu- 



lation of 



Bate 



; 1 2009a I and evolved it to an age of 10 yrs using an 
N-body code under a variety of assumptions regarding the dispersal 
of the molecular cloud. Moeckel & Bate found that the multiplic- 
ity distribution evolved very little during dispersal of the molecular 
cloud and was surprisingly robust to different assumptions regard- 
ing gas dispersal. Even under the assumption of no gas removal at 
lO'^ yrs, although the multiplicities were found to have decreased 
slightly compared with those at the end of the hydrodynamical cal- 
culation, they were still formally consistent. They concluded that 
when star formation occurs in a clustered environment, the multi- 
ple systems that are produced are quite robust against dynamical 
disruption during continued evolution. Therefore, we do not expect 
the multiplicities presented in Fig.[T7]to evolve significantly as the 
stars evolve into a field population. 
In detail, we find: 



Solar-type stars: Duquennoy & Mayor| ^1991^ found an ob- 
served multiplicity fraction of mf = 0.58 ± 0.1. However, the re- 
cent larger survey carried out by |Raghavan et al.]p010| > revised this 
downwards to 0.44 ± 0.02 and they concluded that the higher value 
obtained by [Duquenn oy & Mayor|was due to them overestimating 
their incompleteness correction. The radiation hydrodynamical cal- 
culation gives a multiplicity fraction of 0.42 ± 0.08 over the mass 
range 0.5 — 1.2 Mq which is in good agreement with the result of 
[Raghavan etaLlpOTO) . 

M-dwarfs: [Fischer & Marcyj |T992} found an observed multi- 
plicity fraction of 0.42 ± 0.09. In the mass range 0.1 — 0.5 M© we 
obtain 0.36 ± 0.05. Fischer & Marcy's sample contains stars with 
masses between 0.1 and 0.57 solar masses, but the vast majority 
have masses in the range 0.2 — 0.5 whereas in the simulation 



almost half of the low-mass stars have masses less than 0.2 M©. 
In the 0.2 - 0.5 Mq mass range we obtain 0.38 ± 0.06. All these 
values are consistent with the statistical uncertainties. 



VLM objects: There has been much interest in the multiplicity 
of VLM objects in recent years (Martin et al. 2000 200? 
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gasser et al. 2009, Faherty et al. 2011). For a recent review, see 
|Burgasser et al. ( 2007) . Over the entire mass range of 0.018 — 0.10 
M0, we find a very low multiplicity of just 0.08 ± 0.05, although 
this is twice the value found from the main barotropic calculation 
of [Bate[ ( |2009a[ ). However, the multiplicity drops rapidly with de- 
creasing primary mass and the observed VLM objects tend to have 
high masses. The calculation gives multiplicities of 0.32 ± 0.08 
for the mass range 0.1 — 0.2 Mq, 0.27 ± 0.15 for the mass range 
0.07 - 0.10 Mq, and 0.0 ± 0.05 for the mass range 0.03 - 0.07 
M0. Therefore, to compare with observations it is very important 
to compare like with like. The observed frequency of VLM bina- 
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[Allen| ) (2007| ) estimated the total frequency (including spectroscopic 
systems) to be « 20 — 25%. These surveys typically targeted pri- 
maries with masses in the range 0.03 — 0.1 Mq, but most of these 
objects in fact have masses greater than 0.07 Mq. Thus, the closest 
comparison with our calculation is our frequency of 0.27± 0. 15 for 
the mass range 0.07 — 0.10 Mq. This is in good agreement with ob- 
servations, but the uncertainty is large because of the small number 
of objects. Taking the average over the larger range of 0.03 — 0.20 
Mq gives 0.20 ± 0.05 which is also in good agreement. Because 
the statistics from the radiation hydrodynamical calculation are not 
as good as those of the barotropic simulations it is difficult to deter- 
mine whether or not including radiative feedback has an effect on 
the VLM multiplicity. However, over the range 0.03—0.20 Mq, the 
multiplicity from the rerun barotropic calculation of [Bate[(2009a| ), 
which has the same sink particle accretion radius size as the ra- 
diation hydrodynamical calculation, is 0.17 ± 0.03 which is in 
good agreement with the value obtained from the radiation hydro- 
dynamical calculation. Therefore, the use of radiation hydrodynam- 
ics rather than a barotropic equation of state does not seem to alter 
the VLM multiplicity significantly, and both are in good agreement 
with observations (if small sink particle accretion radii are used; 
[Bate|2009a^ . 

Low-mass brown dwarfs: The frequency of low-mass binary 
brown dwarfs (primary masses less than 30 Jupiter masses) is ob- 
servationally unconstrained. [Bate]j2009a] ) predicted that the multi- 
plicity would continue to fall as the primary mass is decreased and 
that the binary frequency in the mass range 0.01 — 0.03 Mq should 
be ^ 7%. In the radiation hydrodynamical calculation, out of 27 
systems with primary masses < 0.07 Mq there are no multiple 
systems. Thus, although the statistics are not as good, the radiation 
hydrodynamical calculations also predict a very low multiplicity 
for low-mass brown dwarfs. 
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3. 6. 1 Star- VLM binaries 

We turn now to the issue of VLM/brown dwarf companions to stars. 
As in the previous section, we do not consider brown dwarf com- 
panions as such, rather we consider VLM companions (< 0.1 Mq) 
to stars (Jj 0.1 Mq). The radiation hydrodynamical calculation 
produced 8 stellar- VLM systems out of 86 stellar systems, a fre- 
quency of 9 ± 2%. This is indistinguishable from the frequency 
given by the main barotropic calculation (9.0 ± 1.6%). One of the 
8 systems is a 0. 15-M0 star with two VLM companions with with 
semi-major axes of 14 and 194 AU (Table [3|. Another is a 0.52- 
M0 star with a VLM companion at 8.4 AU and a O.19-M0 star at 
20 AU. The other six are binaries. Five of the binaries have primary 
masses in the range 0.1 — 0.2 Mq (semi-major axes of 9, 15, 356, 
721, and 8366 AU), and the sixth is a 0.88 M© star with a 36 M,j 
companion at 620 AU. 

Although the statistics are not as good as the main barotropic 
calculation (which produced 26 stellar- VLM systems), the proper- 
ties are in good agreement. In the main barotropic calculation, 14 
of the primaries had masses between 0.1 — 0.2 Mq, 7 had primary 
masses in the range 0.2 — 0.5 Mq, and 3 had primary masses be- 
tween 0.5 and 0.8 Mq jBate|2009a^ . However, within the statistical 
uncertainties, the frequency of stellar- VLM systems was not found 
to vary with primary mass. Bate (2009a) found a dependence of 
the separations of stellar- VLM binaries on primary mass. He found 
a wide range of separations for primary masses < 0.2 M0, from 
close (< 20 AU) to wide (> 1000 AU) systems, but that the median 
separation increasing strongly with primary mass (< 30 AU for 
stellar masses < 0.2 Mq to ^ 1000 AU for solar- type primaries). 
Although we cannot confirm this trend of median separation with 
primary mass from the radiation hydrodynamical simulation due to 
the smaller numbers of systems, the results are consistent with the 
barotropic results. 

There has been much discussion over the past decade of the 
observed "brown dwarf desert" for close brown dwarf companions 
solar-type stars (frequency « l%; |Marcy & Butler|2000[|Grether &| 
|Lineweaver|2006^ and how this changes for wider separations and 
different primary masses. [McCarthy & Zuckerman| ( |2004| ) found 
that the frequency of wide brown dwarfs to G, K, and M stars be- 
tween 75-300 AU was 1%±1%. The frequencies of wide brown 
dwarf companions to A and B stars ( [Kouwenhoven et al.|2007| ), M 
dwarfs l |Gizis et^L||20031 l, and other brown dwarfs appears to be 
similarly low, although the frequency of wide binary brown dwarfs 
may be higher when they are very young ( [Close et al.|2007^ . Our 
results are consistent with these observations in the sense that we 
do not find brown dwarf companions to solar-type stars in close 
orbits and we find a low frequency of VLM companions at larger 
separations. However, with our small numbers of objects we are 
unable to place strong limits. 

3.6.2 The frequencies of triple and quadruple systems 

Consulting Table |2] we find that the radiation hydrodynamical 
calculation produced 84 single stars/brown dwarfs, 28 binaries, 5 
triples and 7 quadruples. This gives an overall frequency of triple 
and quadruple systems of only 4 ± 2% and 6 ± 2%, respectively. 
These are upper limits because some of these systems may be dis- 
rupted if the calculation were followed longer. Bate (2009a) found 
slightly lower values from the barotropic calculations, although the 
radiation hydrodynamic results agree within the uncertainties. 

For the barotropic calculations, the frequencies of high-order 
multiples were found to increase strongly with primary mass. 



Although the statistics are much poorer for the radiation hy- 
drodynamical calculation, this also appears to be the case. For 
VLM primaries, there are no triples or quadruples out of 38 sys- 
tems. For M-dwarf primaries (0.10 — 0.50 Mq) the frequency of 
triples/quadruples is 8 ± 4%, while for solar-type and intermediate 
mass stars the frequency is « 26 ± 10%. 

Comparing these results with observations, [Fischer & Marcy| 
(1992) found 7 triples and 1 quadruple amongst 99 M-star pri- 
maries giving a frequency of 8 ± 3%, in good agreement. [Raghavan[ 
[et al.[(20T0| found the frequency of triple and higher-order multiple 
systems with solar- type primaries to be 11 ± 1%. For primaries in 
the mass range 0.5 — 1.2 Mq, the radiation hydrodynamical simu- 
lation gives a frequency of 17 ± 12%. In summary, the frequencies 
of triples/quadruples obtained from the radiation hydrodynamical 
calculation are consistent with current observational surveys, but 
the statistical uncertainties are large. 

3.7 Separation distributions of multiples 

The main barotropic calculation of [Bate[ ( [2009a[ ) produced the first 
reasonably large sample of multiple systems from a single hydro- 
dynamical calculation: 58 stellar and 32 VLM binaries, in addition 
to 19 stellar and 4 VLM triple systems and 23 stellar and 2 VLM 
quadruple systems. The radiation hydrodynamical calculation pro- 
duces fewer multiple systems due to the effects of radiative feed- 
back and the because we are not able to follow the evolution as 
far (Table [TJ. However, because the characteristic stellar mass in- 
creases with radiative feedback, it still provides nearly half as many 
multiple stellar systems (25 binaries, 5 triples and 7 quadruples). 
The main difficulty is because the number of VLM objects is more 
than an order of magnitude lower, so the number of VLM multiples 
is very small. The calculation only produces 3 VLM binaries. De- 
spite this, it is of interest to examine the distributions of semi-major 
axes. 

Observationally, the median separation of binaries is found 
to depend on primary mass. [Duquennoy & Mayor[ ( [1991^ found 
that the median separation of solar-type binaries was « 30 AU. In 
the recent larger survey of solar-type stars, [Raghavan et al.lpOlO] ) 
found ~ 40 AU. [Fischer & Marcy| jI992| ) found indications of a 
smaller median separation of « 10 AU for M-dwarf binaries. Fi- 
nally, VLM binaries are found to have a median separation ^ 4 AU 
( [Close et al.|2003| [20071 [Siegler et al.|2005] l, with few VLM bina- 
ries found to have separations greater than 20 AU, particularly in 
the field ( Allen et al. 2007 ). A list of VLM multiple systems can be 
found at http://vlmbinaries.org/ |Close et al.[p007T ) estimated that 
young VLM objects have a wide (> 100 AU) binary frequency of 
~ 6%±3% for ages less than 10 Myr, but only 0.3%±0.I% for 
field VLM objects. 

Although we are able to follow binaries as close as 0.01 Rq 
before they are assumed to merge in the radiation hydrodynami- 
cal calculation, the sink particle accretion radii are 0.5 AU. Thus, 
dissipative interactions between stars and gas are omitted on these 
scales which likely affects the formation of close systems ( [Bate[ 
[et al.|2002a^ . 

In Figure[T8] we present the separation (semi-major axis) dis- 
tributions of the stellar (primary masses greater than 0. 10 Mq) mul- 
tiples. We do not plot the distribution of VLM binaries because 
there are only three systems. The distribution is compared with 
the log-normal distribution from the survey of solar-type stars of 
Raghavan et al. ( [20I0[ ) (which is very similar to that of |Duquennoy[ 
& Mayor|I99I '. The filled histogram gives the separations of bi- 
nary systems, while the double hatched region adds the separations 
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Figure 18. The distributions of separations (semi-major axes) of multiple 
systems with stellar primaries (M* > 0.1 Mq) produced by the radiation 
hydrodynamical calculation. The solid, double hatched, and single hatched 
histograms give the orbital separations of binaries, triples, and quadruples, 
respectively (each triple contributes two separations, each quadruple con- 
tributes three separations). The curve gives the solar-type separation distri- 
bution (scaled to match the area) from Raghavan et al. (2010). The vertical 
dotted line gives the resolution limit of the calculations as determined by 
the accretion radii of the sink particles. 



from triple systems (two separations for each triple, determined 
by decomposing a triple into a binary with a wider companion), 
and the single hatched region includes the separations of quadru- 
ple systems (three separations for each quadruple which may be 
comprised of two binary components or a triple with a wider com- 
panion). The vertical dotted line denotes the sink particle accretion 
radius. 

The median separation (including separations from binary, 
triple, and quadruple systems) of the stellar systems is 15 AU with 
a standard deviation of 0.97 dex (i.e. 1 order of magnitude). Given 
the smaller number of systems, this is in reasonable agreement with 
the value of 26 AU with a standard deviation of 1.15 dex obtained 
by [Batelp009a) for the main barotropic calculation. Both the me- 
dian separation (40 AU) and the dispersion (1.52 dex) obtained by 
[Raghavan et al.| ( [2010^ for solar-type stars are larger. However, for 
the median, it is important to note that most of the primaries from 
the calculation are M-dwarfs, not solar-type stars and )Fischer &| 
|Marcy| ( [l992| l found that M-dwarf binaries have smaller separations 
(median ~ 10 AU). For the width of the distribution, the number 
of close systems is likely underestimated because of the lack of 
dissipation on small scales (see |Bate|2009a^ . The number of wide 
systems is low because the stellar cluster is dense and wide binaries 
can not exist within the cluster There appears to be a similar deficit 
of wide binaries in the Orion Nebula Cluster (Bate, Clarke & Mc-' 
Caughrean 19981 [Scally, Clarke & McCaughrean 1999, Reipurth 



et al.||2007| . Furthermore, [Moeckel & Bate| ( |2010^ and |Kouwen-| 



hoven et S^ ( |2010[ > have shown that wide systems can be formed as 



a star cluster disperses (see also |Moeckel & Clarke|20lT} . 

The three VLM binaries have semi-major axes of 10.6, 26.1, 



and 36.4 AU. With only three systems, no firm conclusions can 
be drawn. However, these semi-major axes are consistent with the 
fact that most observed VLM binaries have projected separations 
^ 20 AU and that wide systems (> 100 AU) are rare. Furthermore, 
|Bate| j2009a| l, who obtained 32 VLM multiples from the main 
barotropic calculation, found that the median separation of VLM 
multiples decreased as the calculation was evolved from « 30 AU 
at 1.04 iff to ~ 10 AU at 1.50 ts because many were still accret- 
ing gas and interacting with other systems early on. He concluded 
that VLM binaries may form with reasonably wide separations and 
evolve to smaller separations (c.f. [Bate et ar]|2002b] l. It is inter- 
esting to note that of the three VLM binaries found here, only the 
closest (10.6 AU) has stopped accreting. Soon after the binary was 
formed at 1.034 tg, a third object formed nearby making it a VLM 
triple. The binary's initial separation was 12.6 AU which grew to 
16.5 AU while the binary was accreting, and then was reduced to 
10.8 AU in a dynamical encounter aXt = 1.10 tn that terminated 
its accretion. The triple survived until another dynamical encounter 
at t = 1 . 17 tff which striped off the wider VLM object and reduced 
the separation of the binary slightly to 10.6 AU. For the other two 
VLM binaries which were still accreting at the end of the simu- 
lation, the 36.4 AU VLM binary was formed at 1.194 tg with a 
separation of 65 AU which decreased continually until the simu- 
lation was stopped. The 26.1 AU binary formed at 1.13 with an 
initial separation of 48.6 AU which quickly decreased to 19.7 AU 
and then grew under the action of accretion to its final value. In 
addition to the theoretical evidence f rom |B ate] j2009a| l that the sep- 
aration distribution of VLM binaries may evolve with time, the ob- 
servational studies oP Close et al.H20d7] l and |Burgasser et aL]j2007| l 
suggest that young wide VLM binaries are disrupted, leading to the 
observed paucity of old wide VLM systems. 

In summary, the radiation hydrodynamical simulation pro- 
duces a stellar separation distribution that is broad with a median 
separation that is in reasonable agreement with field systems. It 
lacks very close systems (presumably due to the lack of dissipa- 
tion on small scales). It also lacks very wide systems, which may 
be formed as the cluster disperses. The VLM binaries are consis- 
tent with the observation that most VLM binaries are close, but with 
only three systems, two of which are still evolving, no stronger con- 
clusions can be drawn. 



3.8 Mass ratios distributions of multiples 

Along with the separation distributions of the multiple systems we 
can investigate the mass ratio distributions. We begin by consider- 
ing only binaries, but we include binaries that are components of 
triple and quadruple systems. A triple system composed of a bi- 
nary with a wider companion contributes the mass ratio from the 
binary, as does a quadruple composed of a triple with a wider com- 
panion. A quadruple composed of two binaries orbiting each other 
contributes two mass ratios — one from each of the binaries. 

Observationally, the mass ratio distribution of binaries also is 
found to depend on primary mass. |Duquennoy & Mayor| ^1991| l 
found that the mass ratio distribution of solar-type binaries peaked 
at M2/M1 ~ 0.2. However, the recent survey of Raghavaii| 
|et al.| ( |2010^ overturns this result. [Raghavan et al.| ( [2010[ > found 
a flat mass ratio distribution for solar-type primaries in the range 
M^jMx — 0.2 — 0.95, with a drop-off at lower mass ratios and 
a strong peak at nearly equal masses (so-called twins; [Tokovinin] 
[2000b| l. They find the mass ratios of pairs in higher-order systems 
follow the same distribution. These results are consistent with the 
earlier study of [Halbwachs et al.| ^2003^ who found a bi-modal 
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Figure 19. The mass ratio distributions of binary systems with stellar primaries in the mass ranges M\ > 0.5 Mq (left) and M\ = 0.1 — 0.5 Mq (centre) and 
VLM primaries (right; Mi < 0.1 M0) produced by the radiation hydrodynamical calculation. The solid black lines give the observed mass ratio distributions 
of |Raghavan et al.|(2010) for pairs with solar-type primaries (left), |Fischer & Marcy ( 1 992 1 for M 1 = 0.3 — O.57M0 (centre, solid line) and A/i = 0.2 — 0.57 
Mq (centre, dashed line), and of the known very-low-mass binary systems from the list at http://vlmbinaries.org/ (right). The observed mass ratio distributions 
have been scaled so that the areas under the distributions (M2/M1 = 0.4 — 1.0 only for the centre panel) match those from the simulation results. The VLM 
binaries produced by the simulation are biased towards equal masses when compared with M-dwarf and solar-type binaries. All three of the VLM binaries 
have M2/M1 > 0.6 while for the M-dwarf binaries the fraction is only 63% and for the more massive primaries the fraction is only 50%. 



distribution for spectroscopic binaries with primary masses in the 
mass range 0.6 — 1.9 Mq and periods Ss 10 years with a broad 
peak in the range M2/M1 =0.2 — 0.7 and a peak for equal-mass 
systems. [Mazeh et al.| ( |2003) found a flat mass ratio distribution for 
spectroscopic binaries with primaries in the mass range 0.6 — 0.85 
Mq. [Fischer & Marcy | ( | 1992) also found a flat mass ratio distribu- 
tion in the range M2/M1 = 0.4 — 1.0 for M-dwarf binaries with 
all periods. In the Taurus- Auriga star-forming region, [Kraus et aL] 
( |2011^ report a flat mass ratio distribution for primaries in the range 
0.7 — 2.5 Mq, but for primaries in the mass range 0.25 — 0.7 Mq 
they find a bias toward equal-mass systems. This change becomes 
even more extreme for VLM binaries, which are found to have a 
strong preference for equal-mass systems ( [Close et al.|2003[[Siegler| 
let al.|2005|[Reid et al.|200"6) . 

In Figure [19] we present the mass ratio distributions of the 
stars with masses 0.5 Mq (left panel), M-dwarfs with masses 
0.1 Mq ^ Ml < 0.5 Mq (centre panel), and VLM objects (right 
panel). We compare the M-dwarf mass ratio distribution to that of 
[Fischer & MarcyH1992^ , and the higher mass stars to the mass ratio 
distribution of pairs with solar-type primaries obtained by Ragha-' 
[van et ar] ( [2010) . The VLM mass ratio distribution is compared 
with the listing of VLM multiples at [http://vlmbinaries.org/[ 

We find that the ratio of near-equal mass systems to systems 
with dissimilar masses decreases going from VLM objects to solar- 
type stars in a similar way to the observed mass ratio distributions, 
although the statistical significance is not strong. Specifically, all 
three of the VLM binaries have A/2 /Mi > 0.6 while for primary 
masses 0.1 — 0.5 Mq the fraction is only 63%, and for solar-type 
stars (> 0.5 Mq) the fraction is 50%. The M-dwarf mass ratio dis- 
tribution is consistent with Fischer & Marcy's distribution. There 
are only three VLM binaries, but they are consistent with the ob- 
servation most VLM systems have mass ratios greater than 0.6. For 
solar-type systems, [Raghavan et al.[pblO[ l obtained a generally flat 
mass ratio distribution, again consistent with the results obtained 
here. 

The barotropic calculations of [Bate| ( [2009a l) also gave higher 
proportions of near equal-mass binaries for VLM binaries than 
M-dwarf binaries (with greater statistical significance), but gave 



similar fractions for solar-type and M-dwarf binaries. At the time, 
'Bate (2009a) concluded that the barotropic calculations did not pro- 
duce enough unequal-mass solar-type binaries because [buquennoyl 
[& Mayor] ( [1991^ found that the mass ratio distribution peaked at 
M2/M1 « 0.2. But the new mass ration distribution obtained by 
Raghavan et al.[ ( |2010^ has reduced the discrepancy. 

As with the VLM separation distribution, Bate (2009a) found 
that the VLM binary mass ratio distribution evolved with time, 
becoming more biased towards equal-mass systems. Both the ap- 
parent evolution of VLM binary separations and mass ratios are 
consistent with the evolution of close binaries discussed by [Bate[ 
[et al.[ ( |2002b| . Dynamical exchange interactions between binaries 
and single objects tend to produce more equal-mass components, 
as does accretion of gas from circumbinary discs or the accretion 
of infalling gas with high specific angular momentum. The evo- 
lution seen in the radiation hydrodynamical calculation is consis- 
tent with this: the 10.6 AU VLM binary had a mass ratio of 0.49 
at 1.05 ts which grew to its final value of 0.61 before the binary 
stopped accreting at 1.10 tfi, while the mass ratios of the two VLM 
binaries which were still accreting at the end of the calculation 
(M2/M1 = 0.94, 0.98) were being equalised by this accretion. 



3.8.1 Mass ratio versus separation 

In Figure |20] we plot mass ratios against separation (semi-major 
axis) for the binaries, triples, and quadruples at the end of the main 
calculation. Note that for this figure we include systems that are 
sub-components of higher-order systems, using filled symbols to 
denote pairs that are binaries (circles), or are components of triples 
(triangles) or quadruples (squares). We also include the mass ratios 
of the wide components of triples and quadruples. 

|Bate| ( |2009a| l found a clear relation between mass ratio and 
separation from the barotropic calculations, with closer binaries 
having a preference for equal masses. He obtained median mass ra- 
tios for fomary separations in the ranges 1—10, 10—100, 100—1000 
and 1000 - 10"* AU of M2/Afi = 0.74, 0.57, 0.68, 0.17, respec- 
tively. Including the mass ratios of triples and quadruples (as de- 
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Figure 20. The mass ratios of binaries (filled circles), pairs in triples (filled 
triangles), pairs in quadruples (filled squared), the wide components of 
triples (open triangles), and the widest components of quadruples (open 
squares) as a function of semi-major axis at the end of the radiation hydro- 
dynamical calculation. For the wide components of triples, the mass ratio 
compares the mass of the widest component to the sum of the masses of 
the two closest components (the pair). For quadruples involving a two bi- 
nary components (pairs), the mass ratio is between the two pairs, and for 
quadruples involving a triple, the mass ratio is between the mass of the 
fourth component and the triple. There is a clear relationship between mass 
ratio and sepai'ation with closer binaries having a greater fraction of near 
equal-mass systems. 



fined in the caption of Figure|20[l, thiese median values became 0.74, 
0.41, 0.15, and 0.07, respectively. 

In the radiation hydrodynamical calculation, we find a similar 
trend (but again, with poorer statistical significance). The median 
binary mass ratios in the separation ranges 1 — 10, 10 — 100 and 
100 - 1000 AU are M2/M1 = 0.62, 0.81, 0.22, respectively. In- 
cluding the mass ratios of triples and quadruples, the median values 
become M^jMx = 0.59,0.74,0.39, respectively. There is only 1 
binary with a separation > 1000 AU. Thus, there seems to be a 
trend for systems with separations ^ 100 AU to have more equal 
masses than for wider systems. A trend of more equal-mass bina- 
ries with decreasing separation is expected from the evolution of 
protobinary systems accreting gas from an envelope (Bate 20001. 
Furthermore, dynamical interactions between binaries and single 
stars tend to tighten binaries at the same time as increasing the bi- 
nary mass ratio through exchange interactions. 

In terms of the so-called twins, the radiation hydrodynam- 
ical calculation produced 43 binaries (including those in triple 
and quadruple systems), of which there are 7 twins (pairs with 
mass ratios M2/M1 > 0.95) and all have semi-major axes less 
than 40 AU. Furthermore, two of these twins are components 
of triple systems. This is in good agreement with observations 
that consistently find that closer binaries have a higher fraction 
of twins (Soderhjelm 19971 |Tokovinin|[2000b[ [Halbwachs eTaTI 
|2003| l. [Tbkovinini ^2000b| lf'ound evidence for the frequency of 
twins falling off for orbital periods greater than 40 days, while 
[Halbwachs et 31.^2003^ found that the fraction of near equal-mass 



systems (M2/M1 > 0.8) is always larger for shorter period bina- 
ries than longer period binaries regardless of the dividing value of 
the period (from just a few days up to 10 years). The most recent 
study o f,Raghavan et al. ( 2010 ) found the mass ratio distribution de- 
pends on period, with less than 1/4 of twins having periods longer 
than 200 years (separations ~ 40 AU) and no twins having separa- 
tions greater than periods of 1000 years (separations of ~ 115 AU). 

Finally, we note Raghavan et al.| {20ro[ l find that more than 
half of the pairs with periods less than 100 days (separation ~ 
0.5 AU) are components of triples, suggesting that dynamical in- 
teractions may be important for their formation (see [Bate et aT] 
|2002b). The closest pair formed in the radiation hydrodynamical 
calculation is a 0.4- AU binary. However, of the 15 pairs with semi- 
major axes less than 5 AU, 8 are binaries while 4 are components of 
triples and 3 are components of quadruples. Thus, the radiation hy- 
drodynamical calculation also results in approximately half of the 
closest pairs being members of higher-order systems. 



3.8.2 Mass ratios of triples and quadruples 

For stellar triple and quadruple systems, |Tokovinin| j2008^ reports 
that triples are observed to have a median outer mass ratio of 0.39 
independent of the outer orbital period while quadruples involving 
two binary sub-components have a similar median outer mass ratio 
of ~ 0.45, but there appears to be a dependence on the outer orbital 
period with systems with shorter outer periods having higher mass 
ratios. Of 9 triple systems, we obtain a median mass ratio of 0.55 
(0.59 excluding triples which are members of quadruple systems). 
There are only 3 quadruple systems consisting of two pairs, all with 
outer mass ratios > 0.8 and outer periods 5.4 < logj^Q(Pd) < 5.9 
(measured in days). [Tokovininl ( [2008[ l finds no outer mass ratios 
< 0.6 for orbital periods logj^Q(Pd) < 5.4 in this orbital period 
range, but a wide range of mass ratios for longer periods. Since we 
only have three systems and they all fall near the apparent observed 
step change it is not possible to draw any firm conclusions. 

IBate|( [2009a[ ) found that quadruples composed of a triple and 
a wide fourth component out number quadruples composed of two 
binaries by 2:1 in the main barotopic calculation. Observationally, 
|Tokoviniri|j2000a| ) finds roughly equal numbers of such quadruples, 
and the radiation hydrodynamical calculation produces a ratio of 
4:3, consistent with the observations. 

In summary, there is no detectable change in the mass ratio 
distributions of binary and higher-order multiple systems when go- 
ing from barotropic to radiation hydrodynamical calculations. In 
both cases, the calculations are consistent with observed trends 
such as a preference for equal-mass binaries when moving to lower 
primary masses and a preference for twins to have close separa- 
tions. 



3.9 Orbital eccentricities 

Observationally, there is an upper envelope to binary eccentrici- 
ties at periods less than about one year, and binaries with periods 
less than 12 days are almost exclusively found to have circular or- 
bits due to tidal circularisation (Duquennoy & Mayor|1991[|Halb-| 



Iwachs et a l.|2003] [Raghavan et al.|2010| >. However, the radiation 
hydrodynamical calculation does not allow us to probe such small 
separations due to the absence of dissipation on scales < 0.5 AU. 
Observations also indicate that eccentricities e < 0.1 are rare for 
periods greater than « 100 days (separations 1 AU). |Ragha-| 
[van et alT| ( [2010^ finds no binaries with e < 0.1 and orbital peri- 
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Figure 21. The eccentricities of binaries (filled circles), pairs in triples (filled triangles), pairs in quadruples (filled squai'ed), the wide components of triples 
(open triangles), and the widest components of quadruples (open squares) as a function of orbital period (left) and mass ratio (right) at the end of the radiation 
hydrodynamical calculation. The distributions look reasonable, but the eight binaries/pairs with the shortest periods all have large eccentricities which may be 
due to the absence of gas dynamics inside 0.5 AU of each sink particle. 



ods greater than 100 days, though they do find that the outer or- 
bits of two triples and one quadruple have e < 0.1. |Duquennoy &] 
|Mayor| ( [l991| l and |Raghavan et ari ( |2010| l also find that the upper- 
eccentricity envelope is dominated by components of triple sys- 
tems, possibly due to the action of the Kozai mechanism l |Kozai| 
|1962| l. Finally, |Halbwachs et al.|(2003[ > find that the eccentricities 
of binaries with mass ratios M2/M1 > 0.8 with periods greater 
than « 10 days (the tidal circularisation radius) are lower than for 
more unequal mass ratio systems. 

In the left panel of Fig. |2T| we plot the eccentricities versus 
orbital period for the binaries, triples and quadruples from the ra- 
diation hydrodynamical calculation. The symbols have the same 
meaning as in Fig. |20] Bate] j2009a| l found that when using sink 
particle radii of 5 AU in barotropic calculations, that there was 
an excess of high eccentricity (e > 0.7) binaries with separations 
< 10 AU. This excess disappeared when the simulation was rerun 
with small accretion radii of 0.5 AU. Following the gas to smaller 
scales allowed dissipative interactions between closer multiple sys- 
tems. Indeed, this was part of the reason that accretion radii of only 
0.5 AU were used for this paper In the radiation hydrodynami- 
cal calculation, the 8 shortest period binaries all have eccentricities 
between 0.6 and 0.8. These are also the 8 closest systems, with 
semi-major axes ranging from 0.4-2.0 AU. Therefore, despite the 
small accretion radii, it is likely that their high eccentricities are 
due, at least in part, to the lack of dissipative interactions with the 
gas. We do note, however, that 3 of the 8 systems are also binary 
components of higher-order systems and, therefore, their high ec- 
centricities may also be related to the observed upper eccentricity 
envelope of binary components of higher-order systems ( [Raghavanl 
let al.|2010| l. 

The mean eccentricity of all 59 orbits is e = 0.35 ± 0.04 with 
a standard deviation of 0.27. The median is e = 0.27. The mean ec- 
centricity of pairs (including components of triples and quadruples) 
is e = 0.38 ± 0.04 with a standard deviation of 0.27. The mean ec- 



centricity of the triples and quadruples is e = 0.26 ± 0.06 with a 
standard deviation of 0.22. The mean eccentricity obtained by |Bate| 
( |2009a| l for the rerun barotropic calculation (with accretion radii of 
0.5 AU) was e — 0.45. The median eccentricity from |Raghavan| 
|et al.lpOTO) is about e = 0.4, so there is reasonable agreement. 

However, [Raghavan et al.| ( |2010^ report a flat distribution of 
eccentricities for periods longer than 12 days out to e = 0.6, 
whereas the radiation hydrodynamical calculation produces more 
than twice as many orbits with e < 0.2 compared to the the in- 
tervals 0.2 ^ e < 0.4 and 0.4 ^ e < 0.6. In particular, there 
are six binaries with e < 0.08 (along with the outer orbits of one 
triple and one quadruple), where as observed systems with e < 0.1 
are rare. Examining Fig. it can be seen that two of the 6 bina- 
ries are the components of triple systems, which may be related to 
the finding of Ragh avan et al.| ( |2010| ) that components of higher- 
order multiple systems can have low eccentricities. Furthermore, 
all but one of these six binaries has a mass ratio M2/M1 > 0.8 
(right panel of Fig. |21| whic h is in qualitative agreement with the 
finding of Halbw ach s et al.| ( |2003| l that near-equal mass binaries 
have smaller eccentricities than more unequal mass ratio systems. 
|Bate|p009a^ also found evidence that near-equal mass binaries had 
smaller eccentricities in the barotropic calculations. In the radia- 
tion hydrodynamical calculation, the median eccentricity of bina- 
ries with mass ratios M2/M1 < 0.8 is e = 0.45 (29 orbits) while 
for M2/M1 > 0.8 the median is e = 0.13 (14 orbits). Excluding 
the 8 shortest period systems (since they likely have high eccentric- 
ities due to the absence of dissipation on small scales) the median 
eccentricity of binaries with mass ratios M2 /M\ < 0.8 is e = 0.42 
(23 orbits) while for M2/M1 > 0.8 the median is e = 0.10 (12 
orbits). Thus, we also find evidence for a link between mass ratio 
and eccentricity such that near-equal mass systems have lower ec- 
centricities, as is observed. As possible explanation for this is that 
accretion, which drives binaries towards equal masses ( [Artymow^ 
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Figure 22. The relative inclinations of the two orbital planes for the 9 triple systems produced by the radiation hydrodynamical calculation. Triples that are 
sub-components of quadruples are plotted as squares. We give plots of the relative orbital orientation angle versus the semi-major axis of the third component 
(left) and versus the period ratio of the long and short period orbits (right). Only the widest triple has a large relative orbital angle. Note that the one system 
with period ratio Pl/^S ~ 3 is likely to be dynamically unstable and to undergo further evolution. 



|icz|1983[|Bate| 19971 [Bate & Bonnell|1997[|Bate|2000| l, may also 
provide dissipation whicli damps eccentricity. 

Finally, we note that VLM binaries are observed to have 
a preference for low eccentricities with a median value of 0.34 
( |Dupuy & Liu|20rT| l. The barotropic calculation of |Bate| ( [2009a^ 
with small accretion radii also produced low-eccentricity VLM bi- 
naries (|Bate 2010b I, with those VLM binaries with separations less 
than 30 AU having a mean eccentricity of 0.23. Unfortunately, the 
radiation hydrodynamical calculation only produces three VLM bi- 
naries. Two of the three do have small eccentricities (the 10.6-AU 
binary has an eccentricity of 0.46, the 26-AU binary has an eccen- 
tricity of 0.013, and the 36-AU binary has an eccentricity of 0.06), 
but they are also still accreting so no firm conclusions can be drawn. 

3.10 Relative alignment of orbital planes for triples 

For a hierarchical triple system there are two orbital planes, one 
corresponding to the short-period orbit and one to the long-period 
orbit. Observationally, it is difficult to determine the relative ori- 
entation angle, of the two orbits of a triple system due to the 
number of quantities that must be measured to fully characterise the 
orbits. However, the mean value of <& can be measured simply from 
knowledge of the number of co-rotating and counter-rotating sys- 
tems ( |Worley|1967||Tokov inin 1993;'Sterzik & Tokovinin |2002| l. 
The first studies i Worley 1967 , van Albada 1968 1 of the rela- 



tive orbital orientations of triple systems found a small tendency 
towards alignment of the angular momentum vectors of the or- 
bits. Of 54 systems with known directions of the relative motions, 
39 showed co-revolution and 15 counter-revolution resulting in a 
mean relative inclination angle of {$) ~ 50°. For 10 visual sys- 
tems with known orbits, 5 systems were found to have $ < 90°, 
2 had $ > 90° and 3 were ambiguous . |Fekel| ( | 1 98 1 1 examined 20 
systems with known orbits and periods of less than 100 years (for 
the wide orbit). He found that 1/3 had non-coplanar orbits. Finally, 



|Sterzik & Tokovinin| ( |2002| l performed the most detailed study to 
date. From 135 visual triple systems for which the relative direc- 
tions of the orbital motions are known they found {(j)) — 67° ± 9° 
and this result was also consistent with 22 systems for which the 
orbits were known. They also found a tendency for the mean rel- 
ative orbital angular momentum angle to increase with increasing 
orbital period ratio (i.e. systems with more similar orbital periods 
tend to be more closely aligned). 

The main barotropic calculation of |Bate| ( [2009a^ produced 
40 triple systems (17 of which were sub-components of quadru- 
ple systems), with a mean relative orbital orientation angle of 
(<1?) — 65° ± 6°, in good agreement with the observed value. 

The radiation hydrodynamical calculation only produced nine 
triple systems, four of which are components of quadruple systems. 
The mean relative orbital orientation angle of the all these triple 
systems is {"!>) = 30° ± 13°, which is about IJa lower than the 
observed value. For the five pure triples, ($) — 36° ± 26°. The 
relative angles are plotted in Fig. |22]as functions of semi-major 
axis and period ratio. It can be seen that all have small relative an- 
gles, except the widest system. We conclude that both the observed 
and simulated triple systems have a tendency towards orbital copla- 
narity, but the small number of systems produced by the radiation 
hydrodynamical calculation prohibits us from making a stronger 
statement. 



3.11 Relative alignment of discs and orbits 

Finally, we consider the relative alignment of the spins of the sink 
particles in binary systems. Unfortunately there is not a direct anal- 
ogy with real binary systems in this case because the sink particles 
are larger than stars and yet smaller than a typical disc. The ori- 
entation of the sink particle spin thus represents the orientation of 
the total angular momentum of the star and the inner part of its 
surrounding disc. This distinction is important because during the 
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Figure 23. The relative inclinations of the rotation axes of the sink pailicles (modelling stars and their inner discs) of the binary systems produced by the 
radiation hydrodynamical calculation as functions of the binary's separation (left) and eccentricity (right). We include binaries that are sub-components of 
triples (triangles) and quadruples (squares). All binaries for which the spins are closely aligned have semi-major axes <^ 30 AU. 



formation of an object the angular momentum usually varies with 
time as gas falls on to it from the turbulent cloud. Thus, the orien- 
tation of the sink particle frequently differs substantially from the 
orientation of its resolved disc (if one exists) and, furthermore, the 
orientations of both the sink particles and their discs change with 
time while the object continues to accrete gas. The orientations may 
evolve with time due to gravitational torques fBate et al.|2000[ l. 

Observationally, [Weis| ( [T9^74 ) found a tendency for alignment 
between the stellar equatorial planes and orbital planes among pri- 
maries in F star binaries, but not A star binaries. The orbital sepa- 
rations were mainly in the 10 — 100 AU range. Similarly, |Guthrie| 
l |198 5) found no correlation for 23 A star binaries with separa- 
tions 10-70 AU. More recently, [Hale| \ 1 994.) considered 73 binary 
and multiple systems containing solar-type stars and found evi- 
dence for approximate coplanarity between the orbital plane and 
the stellar equatorial planes for binary systems with separations 
less than ~ 30 AU and apparently uncorrelated stellar rotation 
and orbital axes for wider systems. For higher-order multiple sys- 
tems, however, non-coplanar systems were found to exist for both 
wide and close orbits. Hale found no evidence to support a dif- 
ference dependent on spectral type, eccentricity or age. In terms 
of circumstellar discs, there is evidence for misaligned discs from 



observations of misaligned jets from protostellar objects i Davis, 



Mundt & Eisloeffel 1994|l, inferred jet precession ( |Eisloffel et al 
1996', 'Da vis et al.|1997| l, and direct observations l| Koresko|1998| 
Stapelfeldt et al.||1998^. However, these are not statistically use- 
ful samples. Monin, Menard & Duchene ( 1998), Donar, Jensen & 
|MathieU| ^1999|), ,Jen sen et al. (2004), Wolf, St ecklum & Henning, 
1 200 1^, and Monin, Menard & Peretto ( 2006 1 used polarimetry to 



study the relative disc alignment in T Tauri wide binary and multi- 
ple systems and all found a preference for disc alignment in bina- 
ries. However, [Jensen et al.| ( |2004l l also found that the wide com- 
ponents of triples and quadruples appear to have random orienta- 
tions. For more massive Herbig Ae/Be binaries, [Baines et al.H2006) 
found that the circumprimary disc was preferentially aligned with 



the orbit and the larger study of | Wheelwright et ar] ( [201 l| l also finds 
that the discs are preferentially aligned with the orbit. 

The barotropic calculations of Bate ( 2009a ) produced ambigu- 
ous results. The main barotropic calculation with large accretion 
radii produced a strong tendency for alignment between sink parti- 
cle spins, but the rerun calculation with smaller accretion radii did 
not show any tendency for alignment l jBate|20l"T| l. 

At the end of the radiation hydrodynamical calculation we plot 
the relative spin angles for the 43 binaries (including those that are 
components of triple and quadruple systems) in Fig.|23]as functions 
of semi-major axis and orbital eccentricity. There are no relative 
spin angles greater than 110°, and only 4 of the 43 systems have 
angles greater than 90°, indicating a strong tendency for alignment. 
The mean relative spin angle is 35° ± 5°, while the median angle is 
20°. For the 28 pure binaries, the mean is 40° ± 6° and the median 
is 27°, while for the binaries that are components of higher-order 
systems the mean is 25° ± 8° and the median is 12°, so there is an 
indication that the spins of binaries that are components of higher- 
order multiples may be more aligned. 

Examining the left panel of Fig.[23]it is clear that the tendency 
for alignment depends on separation: all the binaries that have rel- 
ative spin angles less than 30° have separations less than ~ 30 AU 
or orbital periods less than ~ 400 years. Taking binaries with semi- 
major axes less than 30 AU, the mean relative spin angle is 26° ±5°, 
while those with longer periods have a mean of 70° ± 8° . The right 
panel of Fig.|23]indicates that there may be a weak relation between 
the relative spin angle and the eccentricity, with more circular sys- 
tems having a stronger tendency for alignment. Such a relation may 
come about through accretion onto a binary system or gravitational 
torques between the stars and discs (e.g. |Bate et al.|200 0t, either of 
which would tend to align the components of the binary and may 
damp eccentricity. However, we also note that the distribution of 
relative spin angles seems to be independent of the total mass of 
the binary. 

If the spins of the components of close binaries tend to be 
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aligned with one another, one might also expect the spins to be 
aligned with the orbital plane of the binary. Indeed, this is the case, 
though the alignment is not as strong as for the individual spins. 
Taking the 26 binaries with relative spin angles less than 30°, the 
mean spin-orbit angle is 31° ± 3° with a standard deviation of 19°. 
All of these systems have separations less than 30 AU. For the re- 
maining 17 systems for which the spins are only weakly aligned, 
the mean spin-orbit angle is 48° ± 11° and the standard deviation 
is much larger (64°). 

In summary, for binaries with separations ^ 30 AU, the radia- 
tion hydrodynamical calculation gives strong tendencies for align- 
ment between the spins of the components of binaries and for 
coplanarity of the orbital plane and the equatorial planes of the 
components for binaries. These results are in good agreement with 
the observed coplanarity of observed binaries l |Hale|1994) and in 
qualitative agreement with the many observational studies examin- 
ing disc alignment mentioned above. 



4 CONCLUSIONS 

We have presented results from the largest radiation hydrodynam- 
ical simulation of star cluster formation to date that resolves the 
opacity limit for fragmentation. It also resolves protoplanetary 
discs (radii 1 AU), binaries, and multiple systems. The calcu- 
lation uses sink particles to model the stars and brown dwarfs. We 
discuss in some detail (Section [2. 3[ > the problems associated with 
trying to include the luminosity coming from inside a sink particle's 
accretion radius, concluding that attempts made in the literature to 
date most likely overestimate the luminosity. Although we omit the 
luminosity originating from within each sink particle, we use small 
accretion radii of only 0.5 AU and argue that because protostars 
are observed to be under-luminous (the so-called 'luminosity prob- 
lem') the level of radiative feedback included in the simulation pre- 
sented here may be more realistic than if the extra luminosity was 
included. 

The calculation produced 183 stars and brown dwarfs. This 
number of objects is not as large as that produced from the same 
initial conditions using a barotropic equation of state pate|2009a^ 
because of the impact of radiative feedback. However, it is still 
sufficient to allow comparison of the statistical properties of the 
stars, brown dwarfs and multiple systems with the results of obser- 
vational surveys. Bate (2009a) obtained good agreement between 
observations and barotropic simulations for the properties of multi- 
ple stellar systems, but obtained a brown-dwarf dominated IMF. 
Overall, the radiation hydrodynamical calculation displays good 
agreement with a wide range of observed stellar properties with 
no obvious deficiencies. Together, the barotropic and radiation hy- 
drodynamical calculations imply that the main physical processes 
involved in determining the properties multiple stellar systems are 
gravity and gas dynamics (i.e. dissipative A^-body dynamics), while 
obtaining a realistic IMF also requires radiative feedback. We note, 
however, that the star formation rate in the calculations is much 
higher than observed. To solve this problem may require globally 
unbound molecular clouds and/or the inclusion of magnetic fields 
and kinetic feedback. Our detailed conclusions are as follows. 

(i) The calculation produces an IMF that is statistically indis- 
tinguishable from the parameterisation of the observed IMF by 
IChabrier] P005 [ l , and a ratio of brown dwarfs to stars which is also 
in good agreement with observations. The use of a realistic equa- 
tion of state and radiation hydrodynamics rather than a barotropic 
equation of state decreases the number of brown dwarfs formed 



by an order of magnitude, while having less of an impact on the 
number of stars formed. This corrects the over-production of brown 
dwarfs that is obtained when using a barotropic equation of state. 
We find that the median mass and form of the IMF do not evolve 
significantly during the simulation. 

(ii) As in previous, smaller calculations, the IMF originates 
from competition between accretion and dynamical interactions 
which terminate the accretion and sets an object's final mass. Stars 
and brown dwarfs form the same way, with similar accretion rates 
from the molecular cloud, but stars accrete for longer than brown 
dwarfs before undergoing the dynamical interactions that terminate 
their accretion. The higher characteristic stellar mass that is ob- 
tained when radiative feedback is included comes about because 
the typical distance between objects is larger so that the timescale 
between dynamical interactions is longer and, thus, the objects typ- 
ically accretion to higher masses before their accretion is termi- 
nated. 

(iii) We find that stars have a slightly higher velocity dispersion 
than VLM objects, and binaries have a lower velocity dispersion 
than single objects. 

(iv) We examine the potential effect of dynamical interactions 
on protoplanetary disc sizes. We find that more massive stars have 
had closer encounters. It is difficult to directly associate the closest 
encounter with the radii of protostellar discs because many stars 
accrete new discs after suffering a close encounter. This is partic- 
ularly true for the more massive stars. However, for VLM objects, 
dynamical encounters usually occur soon after their formation and 
terminate their accretion so their truncation radii may more closely 
reflect their disc radii. Under this assumption we find that at least 
20% of VLM objects should have disc radii > 40 AU. In lower 
density star-forming environments this fraction may be expected to 
be larger. 

(v) We find that multiplicity strongly increases with primary 
mass. The results are in good agreement with the observed mul- 
tiplicities of G, K, and M dwarfs and VLM objects. For objects 
with primary masses in the range 0.03 — 0.20 M© the multiplicity 
fraction is 0.20 ± 0.05. We predict that the multiplicity continues 
to drop for lower-mass brown dwarfs. We find very low frequencies 
of VLM companions to stars, in agreement with observations. 

(vi) We examine the separation distributions of binaries, triples 
and quadruples. We find a broad separation distribution for stars 
with a median separation of « 15 AU and a standard deviation of 
1 dex. Unfortunately the calculation only produces three VLM bi- 
naries, two of which are still evolving. However, all of them have 
separations less than 70 AU and the VLM binary that has reached 
its final state has a close separation of 1 1 AU in qualitative agree- 
ment with observations. 

(vii) The mass ratio distributions of solar-type and M-dwarf bi- 
naries are roughly flat, consistent with observations. However, the 
VLM binaries have near-equal masses as appears to be the case for 
observed systems. We find that closer binaries tend to have a higher 
proportion of equal mass components in broad agreement with ob- 
served trends. 

(viii) The eccentricity distribution is broad with no obvious de- 
pendence on period. There may be an excess of short-period highly 
eccentric binaries because of the absence of dissipation on small 
scales due to the use of sink particles. There may also be a weak 
link between mass ratio and eccentricity such that 'twins' have 
lower eccentricities, as is observed. 

(ix) We investigate the relative orientation of the orbital planes 
of triple systems. We find a tendency for orbital alignment, in qual- 
itative agreement with observations. 
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(x) Finally, we study the relative orientations of sink particle 
spins (angular momentum vectors) in binaries (analogous to the ro- 
tation axes of stars and their inner discs). We find that binaries with 
separations ^ 40 AU have a strong tendency for spin alignment, 
in good agreement with observations. We also find that binaries in 
which the spins are closely aligned also have a tendency for align- 
ment of the stellar spins with the orbit. 
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